/*******************************************************************************
*
*  McXtrace, x-ray ray-tracing package
*  Copyright(C) 2007 Risoe National Laboratory.
*
* %I
* Written by: Mads Bertelsen and Erik B Knudsen
* Date: 20.08.15
* Origin: ESS DMSC & DTU Physics
*
* Implementation of a sphere geometry - to be filled with a material
*
* %D
* Part of the Union components, a set of components that work together and thus
*  separates geometry and physics within McXtrace.
* The use of this component requires other components to be used.
*
* 1) One specifies a number of processes using process components
* 2) These are gathered into material definitions using Union_make_material
* 3) Geometries are placed using Union_box/cylinder/sphere, assigned a material
* 4) A Union_master component placed after all of the above
*
* Only in step 4 will any simulation happen, and per default all geometries
*  defined before this master, but after the previous will be simulated here.
*
* There is a dedicated manual available for the Union components
*
* The position of this component is the center of the sphere.
*
* It is allowed to overlap components, but it is not allowed to have two
*  parallel planes that coincide. This will crash the code on run time.
*
*
* %P
* INPUT PARAMETERS:
* radius: [m]          Radius of sphere
* material_string: [str]  material name of this volume, defined using Union_make_material
* priority:   [1]  priority of the volume (can not be the same as another volume) A high priority is on top of low.
* p_interact: [1]  probability to interact with this geometry [0-1]
* visualize:  [1]  set to 0 if you wish to hide this geometry in mcdisplay
* number_of_activations: [1]  Number of subsequent Union_master components that will simulate this geometry
* mask_string: []      Comma separated list of geometry names which this geometry should mask
* mask_setting: []     "All" or "Any", should the masked volume be simulated when the ray is in just one mask, or all.
* target_index:   [1]  Focus on a component this many steps further in the component sequence. 
* target_x:   [m] \
* target_y:   [m] - Position of target to focus at
* target_z:   [m] /
* focus_aw:   [deg] horiz. angular dimension of a rectangular area
* focus_ah:   [deg] vert. angular dimension of a rectangular area
* focus_xw:   [m]   horiz. dimension of a rectangular area
* focus_xh:   [m]   vert. dimension of a rectangular area
* focus_r:    [m]   focusing on circle with this radius
* init:              [string]    Deprecated and unused. Accepted so that instruments written for McStas/McXtrace 3.8.7 and earlier, which name the Union_init component here, still compile.
*
* OUTPUT PARAMETERS:
*
* %L
*
* %E
******************************************************************************/

DEFINE COMPONENT Union_sphere
DEFINITION PARAMETERS ()
SETTING PARAMETERS(string material_string=0, priority, radius, visualize=1, int target_index=0, target_x=0, target_y=0, target_z=0, focus_aw=0, focus_ah=0, focus_xw=0, focus_xh=0, focus_r=0, p_interact=0, string mask_string=0, string mask_setting=0,number_of_activations=1, string init="")
OUTPUT PARAMETERS ()

/* X-ray parameters: (x,y,z,kx,ky,kz ,phi, t, Ex,Ey,Ez,p) */

SHARE
%{
  %include "union-lib"

  void
  mcdisplay_sphere_function (struct lines_to_draw* lines_to_draw_output, int index, struct geometry_struct** Geometries, int number_of_volumes) {
    // Function to call in mcdisplay section of the sample component for this volume
    // One can assume that Geometries[index] refers to a geometry as described in this file
    // The 4 lines describing the sphere are aligned to the local frame of the sphere,
    //   it would be nicer to have them alligned with the global frame so that they show up nicely in
    //   pgplotters on mcdisplay.
    // One could get the current global rotation and use this to counteract this effect.

    double radius = Geometries[index]->geometry_parameters.p_sphere_storage->sph_radius;
    Coords center = Geometries[index]->center;

    Coords direction1 = coords_set (0, 0, 1.0);
    Coords direction2 = coords_set (0, 1.0, 0);
    Coords direction3 = coords_set (1.0, 0, 0);

    struct lines_to_draw lines_to_draw_temp;
    lines_to_draw_temp.number_of_lines = 0;

    lines_to_draw_temp = draw_circle_with_highest_priority (center, direction1, radius, index, Geometries, number_of_volumes, 2);
    merge_lines_to_draw (lines_to_draw_output, &lines_to_draw_temp);

    lines_to_draw_temp = draw_circle_with_highest_priority (center, direction2, radius, index, Geometries, number_of_volumes, 2);
    merge_lines_to_draw (lines_to_draw_output, &lines_to_draw_temp);

    lines_to_draw_temp = draw_circle_with_highest_priority (center, direction3, radius, index, Geometries, number_of_volumes, 2);
    merge_lines_to_draw (lines_to_draw_output, &lines_to_draw_temp);
  };

  void
  initialize_sphere_geometry_from_main_component (struct geometry_struct* sphere) {
    // Function to be called in initialize of the main component
    // This is done as the rotation matrix needs to be relative to the main component instead of global
    // Everything done in initialize in this component file has the rotation matrix relative to global

    // Nothing needs to be done for the sphere
    // If an empty function provides difficulties for some compilers, a dummy operation can be added
    int dummy;
  };

  struct pointer_to_1d_coords_list
  sphere_shell_points (struct geometry_struct* geometry, int max_number_of_points) {
    // Function that returns a number (less than max) of points on the geometry surface
    // If used, remember to free the space allocated.

    int points_per_circle = floor (sqrt (max_number_of_points));
    int number_of_circles = points_per_circle;

    struct pointer_to_1d_coords_list sphere_shell_array;
    sphere_shell_array.elements = malloc (points_per_circle * number_of_circles * sizeof (Coords));
    sphere_shell_array.num_elements = points_per_circle * number_of_circles;

    Coords center = geometry->center;
    double radius = geometry->geometry_parameters.p_sphere_storage->sph_radius;

    Coords direction = coords_set (0, 0, 1.0);

    Rotation rot_matrix;
    rot_set_rotation (rot_matrix, (double)PI / number_of_circles, 0, 0);

    // print_rotation(rot_matrix,"rot matrix");

    // Do the first iteration before the loop to avoid one unecessary matrix operation
    points_on_circle (sphere_shell_array.elements, center, direction, radius, points_per_circle);

    int iterate;
    for (iterate = 1; iterate < number_of_circles; iterate++) {
      direction = rot_apply (rot_matrix, direction);
      points_on_circle (sphere_shell_array.elements + points_per_circle * iterate, center, direction, radius, points_per_circle);
    }
    // Other parts of the program change behavior when I use the above code, may be writing to unwanted parts of memory!

    /*
    // Debug
    for (iterate=0;iterate<sphere_shell_array.num_elements;iterate++) {
      //print_position(sphere_shell_array.elements[iterate],"Sphere shell points");
      //printf("\n%f,%f,%f",sphere_shell_array.elements[iterate].x,sphere_shell_array.elements[iterate].y,sphere_shell_array.elements[iterate].z);
    }
    */

    return sphere_shell_array;
  }

  #ifndef ANY_GEOMETRY_DETECTOR_DECLARE
  #define ANY_GEOMETRY_DETECTOR_DECLARE dummy
  // struct pointer_to_global_geometry_list global_geometry_list = {0,NULL};
  #endif
%}

DECLARE
%{
  // Needed for transport to the main component
  struct global_geometry_element_struct global_geometry_element;

  int loop_index;

  struct Volume_struct this_sphere_volume;
  struct sphere_storage this_sphere_storage;
%}

INITIALIZE
%{
  struct union_state_struct* union_state_p = union_acquire ();
  // Initializes the focusing system for this volume including input sanitation.
  focus_initialize (&this_sphere_volume.geometry, POS_A_COMP_INDEX (INDEX_CURRENT_COMP + target_index), POS_A_CURRENT_COMP, ROT_A_CURRENT_COMP, target_index,
                    target_x, target_y, target_z, focus_aw, focus_ah, focus_xw, focus_xh, focus_r, NAME_CURRENT_COMP);

  // Input sanitation for this geometry
  if (radius <= 0) {
    printf ("\nERROR in Union_sphere named %s, the radius is <= 0. \n", NAME_CURRENT_COMP);
    exit (1);
  }

  // Get global variable from the shared Union state
  struct pointer_to_global_material_list* global_material_list = &union_state_p->u_material_list;
  // Use sanitation
  #ifdef MATERIAL_DETECTOR
  if (global_material_list->num_elements == 0) {
    // Here if the user have defined a material, but only after this material
    printf ("\nERROR: Need to define a material using Union_make_material before using a Union geometry component. \n");
    printf ("       %s was defined before first use of Union_make_material.\n", NAME_CURRENT_COMP);
    exit (1);
  }
  #endif
  #ifndef MATERIAL_DETECTOR
  printf ("\nERROR: Need to define a material using Union_make_material before using a Union geometry component. \n");
  exit (1);
  #endif

  this_sphere_volume.geometry.is_masked_volume = 0;
  this_sphere_volume.geometry.is_exit_volume = 0;
  this_sphere_volume.geometry.is_mask_volume = 0;

  struct pointer_to_global_geometry_list* global_geometry_list = &union_state_p->u_geometry_list;
  // Read the material input, or if it lacks, use automatic linking.
  if (mask_string && strlen (mask_string) && strcmp (mask_string, "NULL") && strcmp (mask_string, "0")) {
    // A mask volume is used to limit the extend of other volumes, called the masked volumes. These are specified in the mask_string.
    // In order for a ray to enter a masked volume, it needs to be both in the region covered by that volume AND the mask volume.
    // When more than
    this_sphere_volume.geometry.mask_mode = 1; // Default is mask mode is ALL
    if (mask_setting && strlen (mask_setting) && strcmp (mask_setting, "NULL") && strcmp (mask_setting, "0")) {
      if (strcmp (mask_setting, "ALL") == 0 || strcmp (mask_setting, "All") == 0)
        this_sphere_volume.geometry.mask_mode = 1;
      else if (strcmp (mask_setting, "ANY") == 0 || strcmp (mask_setting, "Any") == 0)
        this_sphere_volume.geometry.mask_mode = 2;
      else {
        printf ("The mask_mode of component %s is set to %s, but must be either ALL or ANY.\n", NAME_CURRENT_COMP, mask_setting);
        exit (1);
      }
    }

    int found_geometries = 0;
    for (loop_index = 0; loop_index < global_geometry_list->num_elements; loop_index++) {
      // Add mask list
      if (1 == manual_linking_function (global_geometry_list->elements[loop_index].name, mask_string)) {
        add_element_to_int_list (&this_sphere_volume.geometry.mask_list, global_geometry_list->elements[loop_index].component_index);
        add_element_to_int_list (&global_geometry_list->elements[loop_index].Volume->geometry.masked_by_list, INDEX_CURRENT_COMP);
        global_geometry_list->elements[loop_index].Volume->geometry.is_masked_volume = 1;
        if (this_sphere_volume.geometry.mask_mode == 2)
          global_geometry_list->elements[loop_index].Volume->geometry.mask_mode = 2;
        if (this_sphere_volume.geometry.mask_mode == 1) {
          if (global_geometry_list->elements[loop_index].Volume->geometry.is_masked_volume == 1
              && global_geometry_list->elements[loop_index].Volume->geometry.mask_mode != 2)
            // If more than one mask is added to one volume, the ANY mode overwrites the (default) ALL mode.
            global_geometry_list->elements[loop_index].Volume->geometry.mask_mode = 1;
        }

        found_geometries = 1;
      }
    }
    if (found_geometries == 0) {
      printf ("The mask_string in geometry: %s did not find any of the specified volumes in the mask_string %s \n", NAME_CURRENT_COMP, mask_string);
      exit (1);
    }
    this_sphere_volume.p_physics = malloc (sizeof (struct physics_struct));
    this_sphere_volume.p_physics->is_vacuum = 0;                // Makes this volume a vacuum
    this_sphere_volume.p_physics->number_of_processes = (int)0; // Should not be used.
    this_sphere_volume.p_physics->my_a = 0;                     // Should not be used.
    sprintf (this_sphere_volume.p_physics->name, "Mask");
    this_sphere_volume.geometry.is_mask_volume = 1;

    // Read the material input, or if it lacks, use automatic linking.
  } else if (material_string && strlen (material_string) && strcmp (material_string, "NULL") && strcmp (material_string, "0")) {
    // A geometry string was given, use it to determine which material
    if (0 == strcmp (material_string, "vacuum") || 0 == strcmp (material_string, "Vacuum")) {
      // One could have a global physics struct for vacuum instead of creating one for each
      this_sphere_volume.p_physics = malloc (sizeof (struct physics_struct));
      this_sphere_volume.p_physics->is_vacuum = 1; // Makes this volume a vacuum
      this_sphere_volume.p_physics->number_of_processes = (int)0;
      this_sphere_volume.p_physics->my_a = 0; // Should not be used.
      sprintf (this_sphere_volume.p_physics->name, "Vacuum");
    } else if (0 == strcmp (material_string, "exit") || 0 == strcmp (material_string, "Exit")) {
      // One could have a global physics struct for vacuum instead of creating one for each
      this_sphere_volume.p_physics = malloc (sizeof (struct physics_struct));
      this_sphere_volume.p_physics->is_vacuum = 1; // Makes this volume a vacuum
      this_sphere_volume.p_physics->number_of_processes = (int)0;
      this_sphere_volume.p_physics->my_a = 0; // Should not be used.
      this_sphere_volume.geometry.is_exit_volume = 1;
      sprintf (this_sphere_volume.p_physics->name, "Exit");
    } else {
      for (loop_index = 0; loop_index < global_material_list->num_elements; loop_index++) {
        if (0 == strcmp (material_string, global_material_list->elements[loop_index].name)) {
          this_sphere_volume.p_physics = global_material_list->elements[loop_index].physics;
          break;
        }
        if (loop_index == global_material_list->num_elements - 1) {
          printf ("\n");
          printf ("ERROR: The material string \"%s\" in Union geometry \"%s\" did not match a specified material. \n", material_string, NAME_CURRENT_COMP);
          printf ("       The materials available at this point (need to be defined before the geometry): \n");
          for (loop_index = 0; loop_index < global_material_list->num_elements; loop_index++)
            printf ("         %s\n", global_material_list->elements[loop_index].name);
          printf ("\n");
          printf ("       It is also possible to use one of the defualt materials avaiable: \n");
          printf ("           Vacuum (for a Volume without scattering or absorption)\n");
          printf ("           Exit (for a Volume where the ray exits the component if it enters)\n");
          printf ("           Mask (for a Volume that masks existing volumes specified in the mask_string\n");
          exit (1);
        }
      }
    }
  } else {
    // Automatic linking, simply using the last defined material.
    #ifndef MATERIAL_DETECTOR
    printf ("Need to define a material before the geometry to use automatic linking %s.\n", NAME_CURRENT_COMP);
    exit (1);
    #endif
    this_sphere_volume.p_physics = global_material_list->elements[global_material_list->num_elements - 1].physics;
  }

  sprintf (this_sphere_volume.name, NAME_CURRENT_COMP);
  sprintf (this_sphere_volume.geometry.shape, "sphere");
  this_sphere_volume.geometry.eShape = sphere;
  this_sphere_volume.geometry.priority_value = priority;
  // Currently the coordinates will be in absolute space.
  this_sphere_volume.geometry.center = POS_A_CURRENT_COMP;

  this_sphere_volume.geometry.geometry_p_interact = p_interact;
  this_sphere_storage.sph_radius = radius;
  this_sphere_volume.geometry.visualization_on = visualize;
  this_sphere_volume.geometry.geometry_parameters.p_sphere_storage = &this_sphere_storage;
  this_sphere_volume.geometry.within_function = &r_within_sphere;
  this_sphere_volume.geometry.intersect_function = &sample_sphere_intersect;
  this_sphere_volume.geometry.mcdisplay_function = &mcdisplay_sphere_function;
  this_sphere_volume.geometry.initialize_from_main_function = &initialize_sphere_geometry_from_main_component;
  this_sphere_volume.geometry.shell_points = &sphere_shell_points;
  this_sphere_volume.geometry.process_rot_allocated = 0;
  this_sphere_volume.geometry.copy_geometry_parameters = &allocate_sphere_storage_copy;

  rot_copy (this_sphere_volume.geometry.rotation_matrix, ROT_A_CURRENT_COMP);
  rot_transpose (ROT_A_CURRENT_COMP, this_sphere_volume.geometry.transpose_rotation_matrix);

  // Initialize loggers
  this_sphere_volume.loggers.num_elements = 0;
  this_sphere_volume.abs_loggers.num_elements = 0;

  // packing the information into the global_geometry_element, which is then included in the global_geometry_list.
  sprintf (global_geometry_element.name, NAME_CURRENT_COMP);
  global_geometry_element.activation_counter = number_of_activations;
  global_geometry_element.component_index = INDEX_CURRENT_COMP;
  global_geometry_element.Volume = &this_sphere_volume; // Would be nicer if this m was a pointer, now we have the (small) data two places
  add_element_to_geometry_list (global_geometry_list, global_geometry_element);
%}

TRACE
%{
%}

FINALLY
%{
  union_release (NAME_CURRENT_COMP);
%}

END