/*******************************************************************************
*
*  McStas, neutron ray-tracing package
*  Copyright(C) 2007 Risoe National Laboratory.
*
* %I
* Written by: Mads Bertelsen
* Date: 20.08.15
* Origin: University of Copenhagen
*
* Sphere geometry component for the Union components
*
* %D
* Part of the Union components, a set of components that work together and thus
*  sperates geometry and physics within McStas.
* 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 coincides. This will crash the code on run time.
*
*
* %P
* INPUT PARAMETERS:
* radius:                [m]   Radius of sphere
* material_string:       [string]    material name of this volume, defined using Union_make_material
* priority:              [1]   priotiry 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:           [string]    Comma seperated list of geometry names which this geometry should mask
* mask_setting:          [string]    "All" or "Any", should the masked volume be simulated when the ray is in just one mask, or all.
* target_index:          [1]    Focuses on component a component this many steps further in the component sequence
* target_x:              [m]  X position of target to focus at
* target_y:              [m]  Y position of target to focus at
* target_z:              [m]  Z position of target to focus at
* 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]    name of Union_init component (typically "init", default)
*
* CALCULATED PARAMETERS:
*
* %L
*
* %E
******************************************************************************/

DEFINE COMPONENT Union_sphere

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="init")


/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */

SHARE
%{
#ifndef Union
#error "The Union_init component must be included before this Union_sphere component"
#endif


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;
int loop_2_index;
int material_index;

struct Volume_struct this_sphere_volume;
struct sphere_storage this_sphere_storage;
%}

INITIALIZE
%{
// 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);
}

if (_getcomp_index(init) < 0) {
fprintf(stderr,"Union_sphere:%s: Error identifying Union_init component, %s is not a known component name.\n",
NAME_CURRENT_COMP, init);
exit(-1);
}
// Get global variable from init component
struct pointer_to_global_material_list *global_material_list = COMP_GETPAR3(Union_init, init, global_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 = COMP_GETPAR3(Union_init, init, global_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 exit 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,"%s",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,"%s",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);



// TEST Union sphere has been initialized, test shell point function by writing to file.
/*
struct pointer_to_1d_coords_list shell_points;
shell_points = this_sphere_volume.geometry.shell_points(&this_sphere_volume.geometry,400); // using 500 rings with 500 points


FILE *fp;

char filename[124] = "";

strcat(filename,NAME_CURRENT_COMP);
strcat(filename,"_debug.dat");

fp = fopen(filename,"w");

for (loop_index=0;loop_index<shell_points.num_elements;loop_index++)
   fprintf(fp,"%lf %lf %lf\n",shell_points.elements[loop_index].x,shell_points.elements[loop_index].y,shell_points.elements[loop_index].z);

free(shell_points.elements);

fclose(fp);
*/

%}

TRACE
%{
%}

END