/*******************************************************************************
*
*  McStas, neutron ray-tracing package
*  Copyright(C) 2007 Risoe National Laboratory.
*
* %I
* Written by: Mads Bertelsen
* Date: 20.08.15
* Origin: Svanevej 19
*
* A sample component to separate geometry and phsysics
*
* %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 cylinder, and it thus
*  extends yheight/2 up and down along y axis.
*
* 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]          Outer radius volume in (x,z) plane
* yheight: [m]         Cylinder height in (y) direction
* material_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: []      Comma seperated 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.
*
* Focusing options, if left blank, there will be no focusing.
* 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
*
* OUTPUT PARAMETERS:
*
* %L
*
* %E
******************************************************************************/

DEFINE COMPONENT Union_cylinder
DEFINITION PARAMETERS ()
SETTING PARAMETERS(string material_string=0, priority, radius, yheight, 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)
OUTPUT PARAMETERS (loop_index,this_cylinder_volume,global_geometry_element,this_cylinder_storage)

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

SHARE
%{
#ifndef Union
#define Union $Revision: 0.8 $

#include "Union_functions.c"
#include "Union_initialization.c"

#endif


void mcdisplay_cylinder_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 describin the cylinders side are aligned to the local frame of the cylinder,
    //   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 height = Geometries[index]->geometry_parameters.p_cylinder_storage->height;
    double radius = Geometries[index]->geometry_parameters.p_cylinder_storage->cyl_radius;
    Coords direction = Geometries[index]->geometry_parameters.p_cylinder_storage->direction_vector;
    Coords center = Geometries[index]->center;
    
    Coords bottom_point = coords_add(center,coords_scalar_mult(direction,0.5*height));
    Coords top_point = coords_add(center,coords_scalar_mult(direction,-0.5*height));
    
    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(top_point,direction,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(bottom_point,direction,radius,index,Geometries,number_of_volumes,2);
    merge_lines_to_draw(lines_to_draw_output,&lines_to_draw_temp);
    
    Coords point1,point2;
    int iterate,number_of_points=4;
    
    for (iterate=0;iterate<number_of_points;iterate++) {
        point1 = point_on_circle(top_point,direction,radius,iterate,number_of_points);
        point2 = point_on_circle(bottom_point,direction,radius,iterate,number_of_points);
        lines_to_draw_temp = draw_line_with_highest_priority(point1,point2,index,Geometries,number_of_volumes,2);
        merge_lines_to_draw(lines_to_draw_output,&lines_to_draw_temp);
    }
};

void initialize_cylinder_geometry_from_main_component(struct geometry_struct *cylinder) {
    // 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
    
    Coords simple_vector;
    Coords cyl_vector;
    
    // Start with vector that points along the cylinder in the local frame
    simple_vector = coords_set(0,1,0);

    // Rotate the direction vector of the cylinder to the master component frame of reference
    cyl_vector = rot_apply(cylinder->rotation_matrix,simple_vector);
    NORM(cyl_vector.x,cyl_vector.y,cyl_vector.z);
    cylinder->geometry_parameters.p_cylinder_storage->direction_vector.x = cyl_vector.x;
    cylinder->geometry_parameters.p_cylinder_storage->direction_vector.y = cyl_vector.y;
    cylinder->geometry_parameters.p_cylinder_storage->direction_vector.z = cyl_vector.z;
    // if (verbal == 1) printf("Cords vector1 = (%f,%f,%f)\n",cyl_vector.x,cyl_vector.y,
}

struct pointer_to_1d_coords_list cylinder_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(max_number_of_points/2.0);
  
  struct pointer_to_1d_coords_list cylinder_shell_array;
  cylinder_shell_array.elements = malloc(2*points_per_circle*sizeof(Coords));
  cylinder_shell_array.num_elements = 2*points_per_circle;
  
  Coords cyl_direction = geometry->geometry_parameters.p_cylinder_storage->direction_vector;
  Coords center = geometry->center;
  double radius = geometry->geometry_parameters.p_cylinder_storage->cyl_radius;
  double height = geometry->geometry_parameters.p_cylinder_storage->height;
  
  Coords cyl_top_point = coords_add(center,coords_scalar_mult(cyl_direction,0.5*height));
  Coords cyl_bottom_point = coords_add(center,coords_scalar_mult(cyl_direction,-0.5*height));
  
  points_on_circle(cylinder_shell_array.elements,cyl_top_point,cyl_direction,radius,points_per_circle);
  // Need to verify this pointer arithimatic works as intended
  points_on_circle(cylinder_shell_array.elements+points_per_circle,cyl_bottom_point,cyl_direction,radius,points_per_circle);
  
  return cylinder_shell_array;
}

%}

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

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

int dummy;
int loop_index,found_geometries;
int loop_2_index;
int material_index;

struct Volume_struct this_cylinder_volume;
struct cylinder_storage this_cylinder_storage;
%}

INITIALIZE
%{
// Initializes the focusing system for this volume including input sanitation.
focus_initialize(&this_cylinder_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_cylinder named %s, the radius is <= 0. \n",NAME_CURRENT_COMP);
  exit(1);
}

if (yheight <= 0) {
  printf("\nERROR in Union_cylinder named %s, yheight is <= 0. \n",NAME_CURRENT_COMP);
  exit(1);
}

// 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_cylinder_volume.geometry.is_masked_volume = 0;
this_cylinder_volume.geometry.is_exit_volume = 0;
this_cylinder_volume.geometry.is_mask_volume = 0;

// 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_cylinder_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_cylinder_volume.geometry.mask_mode = 1;
        else if (strcmp(mask_setting,"ANY") == 0 || strcmp(mask_setting,"Any") == 0) this_cylinder_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);
        }
    }
    
    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_cylinder_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_cylinder_volume.geometry.mask_mode == 2)
                global_geometry_list.elements[loop_index].Volume->geometry.mask_mode = 2;
            if (this_cylinder_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_cylinder_volume.p_physics = malloc(sizeof(struct physics_struct));
    this_cylinder_volume.p_physics->is_vacuum = 0; // Makes this volume a vacuum
    this_cylinder_volume.p_physics->number_of_processes = (int) 0; // Should not be used.
    this_cylinder_volume.p_physics->my_a = 0; // Should not be used.
    sprintf(this_cylinder_volume.p_physics->name,"Mask");
    this_cylinder_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_cylinder_volume.p_physics = malloc(sizeof(struct physics_struct));
        this_cylinder_volume.p_physics->is_vacuum = 1; // Makes this volume a vacuum
        this_cylinder_volume.p_physics->number_of_processes = (int) 0;
        this_cylinder_volume.p_physics->my_a = 0; // Should not be used.
        sprintf(this_cylinder_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_cylinder_volume.p_physics = malloc(sizeof(struct physics_struct));
        this_cylinder_volume.p_physics->is_vacuum = 1; // Makes this volume a vacuum
        this_cylinder_volume.p_physics->number_of_processes = (int) 0;
        this_cylinder_volume.p_physics->my_a = 0; // Should not be used.
        this_cylinder_volume.geometry.is_exit_volume = 1;
        sprintf(this_cylinder_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_cylinder_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_cylinder_volume.p_physics = global_material_list.elements[global_material_list.num_elements-1].physics;
}

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

this_cylinder_volume.geometry.geometry_p_interact = p_interact;
this_cylinder_storage.cyl_radius = radius;
this_cylinder_storage.height = yheight;
this_cylinder_volume.geometry.visualization_on = visualize;
this_cylinder_volume.geometry.geometry_parameters.p_cylinder_storage = &this_cylinder_storage;
this_cylinder_volume.geometry.within_function = &r_within_cylinder;
this_cylinder_volume.geometry.intersect_function = &sample_cylinder_intersect;
this_cylinder_volume.geometry.mcdisplay_function = &mcdisplay_cylinder_function;
this_cylinder_volume.geometry.shell_points = &cylinder_shell_points;
this_cylinder_volume.geometry.initialize_from_main_function = &initialize_cylinder_geometry_from_main_component;
this_cylinder_volume.geometry.process_rot_allocated = 0;
this_cylinder_volume.geometry.copy_geometry_parameters = &allocate_cylinder_storage_copy;

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

// Initialize loggers
this_cylinder_volume.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_cylinder_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
%{
dummy = 1;
%}

END