/*******************************************************************************
*
*  McStas, neutron ray-tracing package
*  Copyright(C) 2007 Risoe National Laboratory.
*
* %I
* Written by: Martin Olsen
* On a template by: Mads Bertelsen
* Date: 17.09.18
* Origin: University of Copenhagen
*
* Component for including 3D mesh in Union geometry
*
* %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 mesh component loads a 3D stl file as the geometry. The mesh geometry
* can not yet be overlapped with the basic geometries (sphere, cylinder, box, cone),
* but can be overlapped with other mesh geometries.
*
* 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:
* filename: [str] Name of stl file that contains the 3D geometry
* 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]
* 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] name of Union_init component (typically "init", default)
*
* CALCULATED PARAMETERS:
*
* %L
*
* %E
******************************************************************************/

DEFINE COMPONENT Union_mesh

SETTING PARAMETERS(string filename = 0,string material_string=0, priority, 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_mesh component"
#endif



Coords get_coords_from_string(char *line){

int i;

// Find index first number
int first_number_index = 0;
int chars_in_line = strlen(line);
//printf("\n length of line = %i",chars_in_line);

for (i=0 ; i < chars_in_line ; i++ ) {
    ////printf("\n line[%i] = %i",i,line[i]);
    if (line[i] > 44 && line[i] < 58 && first_number_index == 0 && line[i-1] == 32) {
        first_number_index = i;
        i = chars_in_line;
    }
}

int space_index[2];
int counter = 0;
// Find index of spaces
for (i=first_number_index ; i < chars_in_line ; i++ ) {
    ////printf("\n line[%i] = %i",i,line[i]);
    if (line[i] == 32) {
        space_index[counter] = i;
        ////printf("\n space_index[%i] = %i",i,space_index[counter]);
        counter ++;
    }

}

char dest[15];

// Assign x
////printf("\n %s ",line+(first_number_index));
////printf("\n %i ",(space_index[0]-first_number_index));
//strncpy(dest, line+(first_number_index), (space_index[0]-first_number_index));
for (i = 0 ; i < space_index[0]-first_number_index ; i++ ){
    dest[i] = line[first_number_index + i];
}
////printf("\n %s ",dest);
double X = atof(dest);

for (i = 0 ; i < space_index[1]-space_index[0] ; i++ ){
    dest[i] = line[space_index[0] + 1 + i];
}
////printf("\n %s ",dest);
double Y = atof(dest);


for (i = 0 ; i < chars_in_line-space_index[1] ; i++ ){
    dest[i] = line[space_index[1]  + i];
    //printf("\n Current index is %i, max is %i",space_index[1] + 1 + i,chars_in_line);
}
////printf("\n %s ",dest);
double Z = atof(dest);

Coords fromString = coords_set(X,Y,Z);

return fromString;

};

int coord_comp(Coords A,Coords B) {
    if (A.x==B.x && A.y==B.y && A.z==B.z){
        return 1;
    }
    return 0;
};



void mcdisplay_mesh_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

    printf("\nTEST 0");
    int n_facets = Geometries[index]->geometry_parameters.p_mesh_storage->n_facets;
	double *v1_x = Geometries[index]->geometry_parameters.p_mesh_storage->v1_x;
	double *v1_y = Geometries[index]->geometry_parameters.p_mesh_storage->v1_y;
	double *v1_z = Geometries[index]->geometry_parameters.p_mesh_storage->v1_z;
	double *v2_x = Geometries[index]->geometry_parameters.p_mesh_storage->v2_x;
	double *v2_y = Geometries[index]->geometry_parameters.p_mesh_storage->v2_y;
	double *v2_z = Geometries[index]->geometry_parameters.p_mesh_storage->v2_z;
	double *v3_x = Geometries[index]->geometry_parameters.p_mesh_storage->v3_x;
	double *v3_y = Geometries[index]->geometry_parameters.p_mesh_storage->v3_y;
	double *v3_z = Geometries[index]->geometry_parameters.p_mesh_storage->v3_z;

    Coords center = Geometries[index]->center;

    struct lines_to_draw lines_to_draw_temp;
    lines_to_draw_temp.number_of_lines = 0;
    
    Coords point1,point2,point3;
    int iterate, i,j;
    int print1 = 0;
    int print2 = 0;
    int print3 = 0;

    Coords list_startpoints[n_facets*3];
    Coords list_endpoints[n_facets*3];

	int counter=0;
    // For every triangle it should ad three lines    
    for (iterate=0 ; iterate<n_facets ; iterate++) {
        point1 = coords_add(rot_apply(Geometries[index]->rotation_matrix,coords_set(*(v1_x+iterate),*(v1_y+iterate),*(v1_z+iterate))),center);
        point2 = coords_add(rot_apply(Geometries[index]->rotation_matrix,coords_set(*(v2_x+iterate),*(v2_y+iterate),*(v2_z+iterate))),center);
	    point3 = coords_add(rot_apply(Geometries[index]->rotation_matrix,coords_set(*(v3_x+iterate),*(v3_y+iterate),*(v3_z+iterate))),center);

	    print1 = 1;
    	print2 = 1;
    	print3 = 1;
	
	// Make sure it does not print a line if it is already printed.... (might take a while?)
	for (i = 0 ; i < counter ; i++){
		if (print1 == 1 && coord_comp(point1 , list_startpoints[i])){
			for (j = 0 ; j < counter  ; j++){
				if (coord_comp(point2 , list_startpoints[i])){
					print1 = 0;
				}
			}
		}
		if (print2 == 1 && coord_comp(point2 , list_startpoints[i])){
			for (j = 0 ; j < counter  ; j++){
				if (coord_comp(point1 , list_startpoints[i])){
					print1 = 0;
				}
			}
		}
		if (print2 == 1 && coord_comp(point2 , list_startpoints[i]) ){
			for (j = 0 ; j < counter  ; j++){
				if (coord_comp(point3 , list_startpoints[i])){
					print2 = 0;
				}
			}
		}
		if (print3 == 1 && coord_comp(point3 , list_startpoints[i]) ){
			for (j = 0 ; j < counter  ; j++){
				if (coord_comp(point2 , list_startpoints[i])){
					print2 = 0;
				}
			}
		}
		if (print1 == 1 && coord_comp(point1 , list_startpoints[i]) ){
			for (j = 0 ; j < counter ; j++){
				if (coord_comp(point1 , list_startpoints[i])){
					print3 = 0;
				}
			}
		}
		if (print3 == 1 && coord_comp(point3 , list_startpoints[i])){
			for (j = 0 ; j < counter  ; j++){
				if (coord_comp(point1 , list_startpoints[i])){
					print3 = 0;
				}
			}
		}

	}
	
	
	// Create lines	
	// Line 1
	if (print1 == 1){
		lines_to_draw_temp = draw_line_with_highest_priority(point1,point2,index,Geometries,number_of_volumes,100);
		merge_lines_to_draw(lines_to_draw_output,&lines_to_draw_temp);
		list_startpoints[counter] = point1;
		list_endpoints[counter] = point2;
		counter++;
	}
	// Line 2
	if (print2 == 1){
		lines_to_draw_temp = draw_line_with_highest_priority(point2,point3,index,Geometries,number_of_volumes,100);
		merge_lines_to_draw(lines_to_draw_output,&lines_to_draw_temp);
		list_startpoints[counter] = point2;
		list_endpoints[counter] = point3;
		counter++;
	}
	// Line 3
	if (print3 == 1){
		lines_to_draw_temp = draw_line_with_highest_priority(point3,point1,index,Geometries,number_of_volumes,100);
	 	merge_lines_to_draw(lines_to_draw_output,&lines_to_draw_temp);
		list_startpoints[counter] = point3;
		list_endpoints[counter] = point1;
		counter++;
	}
	
    }
};

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

    // Rotate the direction vector of the mesh to the master component frame of reference
    mesh_vector = rot_apply(mesh->rotation_matrix,simple_vector);
    NORM(mesh_vector.x,mesh_vector.y,mesh_vector.z);
    mesh->geometry_parameters.p_mesh_storage->direction_vector.x = mesh_vector.x;
    mesh->geometry_parameters.p_mesh_storage->direction_vector.y = mesh_vector.y;
    mesh->geometry_parameters.p_mesh_storage->direction_vector.z = mesh_vector.z;
    // if (verbal == 1) printf("Cords vector1 = (%f,%f,%f)\n",mesh_vector.x,mesh_vector.y,
    
    mesh->geometry_parameters.p_mesh_storage->Bounding_Box_Center = rot_apply(mesh->rotation_matrix, mesh->geometry_parameters.p_mesh_storage->Bounding_Box_Center);
    
    /*
    // Works for pure translation
    print_position(mesh->geometry_parameters.p_mesh_storage->Bounding_Box_Center, "BB before adjustment");
    mesh->geometry_parameters.p_mesh_storage->Bounding_Box_Center = coords_add(mesh->geometry_parameters.p_mesh_storage->Bounding_Box_Center, mesh->center);
    print_position(mesh->geometry_parameters.p_mesh_storage->Bounding_Box_Center, "BB after adjustment");
    */
}

struct pointer_to_1d_coords_list mesh_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
  // Run trhough all points in list of faces, and remove dublicates
  // There are three points in a face and very often these will be dublicated a few times. This removes dublicates to boost performance down stream...
  
  
  struct pointer_to_1d_coords_list mesh_shell_array;
  
	int n_facets = geometry->geometry_parameters.p_mesh_storage->n_facets;
	double *v1_x = geometry->geometry_parameters.p_mesh_storage->v1_x;
	double *v1_y = geometry->geometry_parameters.p_mesh_storage->v1_y;
	double *v1_z = geometry->geometry_parameters.p_mesh_storage->v1_z;
	double *v2_x = geometry->geometry_parameters.p_mesh_storage->v2_x;
	double *v2_y = geometry->geometry_parameters.p_mesh_storage->v2_y;
	double *v2_z = geometry->geometry_parameters.p_mesh_storage->v2_z;
	double *v3_x = geometry->geometry_parameters.p_mesh_storage->v3_x;
	double *v3_y = geometry->geometry_parameters.p_mesh_storage->v3_y;
	double *v3_z = geometry->geometry_parameters.p_mesh_storage->v3_z;
	int number_of_points_in_array = 0;
	mesh_shell_array.elements = malloc(3*n_facets * sizeof(Coords));
	int is_dublicate = 0;
	Coords this_vert;
	int i,j;
	
	printf("\n CREATE SHELL POINTS");
	printf("\n n_verts (likely dublicated) = %i",n_facets*3);
	for (i=0 ; i < n_facets ; i++){
		
		// v1
		is_dublicate = 0;
		this_vert = coords_set(*(v1_x+i),*(v1_y+i),*(v1_z+i));
		//printf("\n this point [%f,%f,%f] ",this_vert.x,this_vert.y,this_vert.z);
		
		// test if dublicate
		for (j = 0; j < number_of_points_in_array ; j++ ){
			if (this_vert.x == mesh_shell_array.elements[j].x && this_vert.y == mesh_shell_array.elements[j].y && this_vert.z == mesh_shell_array.elements[j].z){
				is_dublicate = 1;
				//printf("\n point [%f,%f,%f] is a dublicate of [%f,%f,%f]",this_vert.x,this_vert.y,this_vert.z,mesh_shell_array.elements[j].x,mesh_shell_array.elements[j].y,mesh_shell_array.elements[j].z);
				j = number_of_points_in_array;
			}
		}
		if (is_dublicate == 0){
			mesh_shell_array.elements[number_of_points_in_array] = this_vert;
			number_of_points_in_array += 1;
			//printf("\n Added a point [%f,%f,%f] ",this_vert.x,this_vert.y,this_vert.z);
		}
	
	
	
		// v2
		is_dublicate = 0;
		this_vert = coords_set(*(v2_x+i),*(v2_y+i),*(v2_z+i));
		
		// test if dublicate
		for (j = 0; j < number_of_points_in_array ; j++){
			if (this_vert.x == mesh_shell_array.elements[j].x && this_vert.y == mesh_shell_array.elements[j].y && this_vert.z == mesh_shell_array.elements[j].z){
				is_dublicate = 1;
				//printf("\n point [%f,%f,%f] is a dublicate of [%f,%f,%f]",this_vert.x,this_vert.y,this_vert.z,mesh_shell_array.elements[j].x,mesh_shell_array.elements[j].y,mesh_shell_array.elements[j].z);
				j = number_of_points_in_array;
			}
		}
		if (is_dublicate == 0){
			mesh_shell_array.elements[number_of_points_in_array] = this_vert;
			number_of_points_in_array += 1;
			//printf("\n Added a point [%f,%f,%f] ",this_vert.x,this_vert.y,this_vert.z);
		}
	
	
	
		// v3
		is_dublicate = 0;
		this_vert = coords_set(*(v3_x+i),*(v3_y+i),*(v3_z+i));
		
		// test if dublicate
		for (j = 0; j < number_of_points_in_array ; j++ ){
			if (this_vert.x == mesh_shell_array.elements[j].x && this_vert.y == mesh_shell_array.elements[j].y && this_vert.z == mesh_shell_array.elements[j].z){
				is_dublicate = 1;
				//printf("\n point [%f,%f,%f] is a dublicate of [%f,%f,%f]",this_vert.x,this_vert.y,this_vert.z,mesh_shell_array.elements[j].x,mesh_shell_array.elements[j].y,mesh_shell_array.elements[j].z);
				j = number_of_points_in_array;
			}
		}
		if (is_dublicate == 0){
			mesh_shell_array.elements[number_of_points_in_array] = this_vert;
			number_of_points_in_array += 1;
			//printf("\n Added a point [%f,%f,%f] ",this_vert.x,this_vert.y,this_vert.z);
		}
	}
    
  printf("\n test is there to high j? j = %i, number_of_points_in_array = %i , lenght of array = %i ",j,number_of_points_in_array,3*n_facets * sizeof(Coords));
  j = number_of_points_in_array - 1; // Last legal index, currently j is out of bounds.
  printf("\n test= [%f,%f,%f]",mesh_shell_array.elements[j].x ,mesh_shell_array.elements[j].y,mesh_shell_array.elements[j].z);
  
  mesh_shell_array.num_elements = number_of_points_in_array;
  //free(mesh_shell_array.elements);
  printf("\n SHELL POINTS: DONE");
  printf("\n SHELL POINTS: created %i shell points in mesh",mesh_shell_array.num_elements);
    //printf("\n test last element = [%f,%f,%f]",mesh_shell_array.elements[mesh_shell_array.num_elements-1].x ,mesh_shell_array.elements[mesh_shell_array.num_elements-1].y,mesh_shell_array.elements[mesh_shell_array.num_elements-1].z);
  return mesh_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_mesh_volume;
struct mesh_storage this_mesh_storage;
%}

INITIALIZE
%{
// Initializes the focusing system for this volume including input sanitation.
focus_initialize(&this_mesh_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);



if (_getcomp_index(init) < 0) {
fprintf(stderr,"Union_mesh:%s: Error identifying Union_init component, %s is not a known component name.\n",
NAME_CURRENT_COMP, init);
exit(-1);
}
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_mesh_volume.geometry.is_masked_volume = 0;
this_mesh_volume.geometry.is_exit_volume = 0;
this_mesh_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_mesh_volume.geometry.mask_mode = 1; // Default 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_mesh_volume.geometry.mask_mode = 1;
        else if (strcmp(mask_setting,"ANY") == 0 || strcmp(mask_setting,"Any") == 0) this_mesh_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_mesh_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_mesh_volume.geometry.mask_mode == 2)
                global_geometry_list->elements[loop_index].Volume->geometry.mask_mode = 2;
            if (this_mesh_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_mesh_volume.p_physics = malloc(sizeof(struct physics_struct));
    this_mesh_volume.p_physics->is_vacuum = 0; // Makes this volume a vacuum
    this_mesh_volume.p_physics->number_of_processes = (int) 0; // Should not be used.
    this_mesh_volume.p_physics->my_a = 0; // Should not be used.
    sprintf(this_mesh_volume.p_physics->name,"Mask");
    this_mesh_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_mesh_volume.p_physics = malloc(sizeof(struct physics_struct));
        this_mesh_volume.p_physics->is_vacuum = 1; // Makes this volume a vacuum
        this_mesh_volume.p_physics->number_of_processes = (int) 0;
        this_mesh_volume.p_physics->my_a = 0; // Should not be used.
        sprintf(this_mesh_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_mesh_volume.p_physics = malloc(sizeof(struct physics_struct));
        this_mesh_volume.p_physics->is_vacuum = 1; // Makes this volume a vacuum
        this_mesh_volume.p_physics->number_of_processes = (int) 0;
        this_mesh_volume.p_physics->my_a = 0; // Should not be used.
        this_mesh_volume.geometry.is_exit_volume = 1;
        sprintf(this_mesh_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_mesh_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_mesh_volume.p_physics = global_material_list->elements[global_material_list->num_elements-1].physics;
}


// Read input file and put into storage 
// test 2
	printf("\n START read file: ");

	//char *filename = "text.stl";
	FILE *fp;
	fp = fopen(filename,"r");
	int n_lines;
	char  buffer[250];
	printf("\n test1");
	n_lines = 0;
	  while (fgets(buffer, sizeof buffer, fp)) {
		++n_lines;
	  }
	  fclose(fp);
	printf("\n n_lines: %i",n_lines);
	printf("\n n_facets: %i",(n_lines-2)/7);
	int n_facets = (n_lines-2)/7;
	
	printf("\n done");
	//printf("\n n_facets: %i",n_facets);
	
	
    char *line = NULL;
    size_t len = 0;
	char dest[9]; // 8 chars + terminator
	int iter = 0; int counter = 0;
    ssize_t read;
	
	Coords normal[n_facets+1];
	Coords v1[n_facets+1];
	Coords v2[n_facets+1];
	Coords v3[n_facets+1];
	
	//#include "MeshFunctions/get_coords_from_string.c"
	printf("\nTEST n_facets = %i\n",n_facets);
	Coords fromString;
    fp = fopen(filename, "r");
    if (fp == NULL)
        exit(EXIT_FAILURE);
    while ((read = getline(&line, &len, fp)) != -1) {
		//printf("______________________-----NEW------_____________________\n");
		//printf("%s\n", line);
		
			if (iter==1){
				if (line[0] == 101){
					//printf("\nDone reading file\n");
				}else{
					
					// Normal vector
					fromString = get_coords_from_string(line);
					normal[counter]=fromString;
				
					//printf("\n normal: [%i]: normal[%i] = [%f,%f,%f] \n",iter,counter,normal[counter].x,normal[counter].y,normal[counter].z);
				}
			}
			if (iter==3){
				fromString = get_coords_from_string(line);
				v1[counter]=fromString;
				//printf("\n v1 [%i]: v1[%i] = [%f,%f,%f] \n",iter,counter,v1[counter].x,v1[counter].y,v1[counter].z);
			}
			if (iter==4){
				fromString = get_coords_from_string(line);
				v2[counter]=fromString;
				//printf("\n v2 [%i]: v2[%i] = [%f,%f,%f] \n",iter,counter,v2[counter].x,v2[counter].y,v2[counter].z);
			}
			if (iter==5){
				fromString = get_coords_from_string(line);
				v3[counter]=fromString;
				//printf("\n v3 [%i]: v3[%i] = [%f,%f,%f] \n",iter,counter,v3[counter].x,v3[counter].y,v3[counter].z);
				
			}
			if (iter==7){
				iter = 0;
				counter ++;
			}
		//printf("\n TESTTESTTEST:  [%f,%f,%f] \n",this_mesh_storage.poligon_list[counter].normal.x,this_mesh_storage.poligon_list[counter].normal.y,this_mesh_storage.poligon_list[counter].normal.z);
		
		iter++; 
    }
	
    fclose(fp);
	
	int i;
	
	// Transform coordinate (translate then rotate)
	printf("\n\nPOS_A_CURRENT_COMP [%f,%f,%f]",POS_A_CURRENT_COMP.x,POS_A_CURRENT_COMP.y,POS_A_CURRENT_COMP.z);
	for (i = 0 ; i < counter ; i++ ){
		// run through existing arrays without creating new ones. This is a very bad habbit!  PLZ fix
		
        /*
		// v1
		v1[i].x = v1[i].x + POS_A_CURRENT_COMP.x;
		v1[i].y = v1[i].y + POS_A_CURRENT_COMP.y;
		v1[i].z = v1[i].z + POS_A_CURRENT_COMP.z;
		//v1[i] = rot_apply(this_mesh_volume.geometry.transpose_rotation_matrix,v1[i]);

		// v2
		v2[i].x = v2[i].x + POS_A_CURRENT_COMP.x;
		v2[i].y = v2[i].y + POS_A_CURRENT_COMP.y;
		v2[i].z = v2[i].z + POS_A_CURRENT_COMP.z;
		//v2[i] = rot_apply(this_mesh_volume.geometry.transpose_rotation_matrix,v2[i]);
		
		// v3
		v3[i].x = v3[i].x + POS_A_CURRENT_COMP.x;
		v3[i].y = v3[i].y + POS_A_CURRENT_COMP.y;
		v3[i].z = v3[i].z + POS_A_CURRENT_COMP.z;
		//v3[i] = rot_apply(this_mesh_volume.geometry.transpose_rotation_matrix,v3[i]);
        */
        
        /*
        // v1
        v1[i].x = v1[i].x + POS_A_CURRENT_COMP.x;
        v1[i].y = v1[i].y + POS_A_CURRENT_COMP.y;
        v1[i].z = v1[i].z + POS_A_CURRENT_COMP.z;
        //v1[i] = rot_apply(this_mesh_volume.geometry.transpose_rotation_matrix,v1[i]);

        // v2
        v2[i].x = v2[i].x + POS_A_CURRENT_COMP.x;
        v2[i].y = v2[i].y + POS_A_CURRENT_COMP.y;
        v2[i].z = v2[i].z + POS_A_CURRENT_COMP.z;
        //v2[i] = rot_apply(this_mesh_volume.geometry.transpose_rotation_matrix,v2[i]);
        
        // v3
        v3[i].x = v3[i].x + POS_A_CURRENT_COMP.x;
        v3[i].y = v3[i].y + POS_A_CURRENT_COMP.y;
        v3[i].z = v3[i].z + POS_A_CURRENT_COMP.z;
        */
		
	}
	printf("\n TEST inf v1 is long enough! i = %i , n_facets = %i, counter = %i \n",i,n_facets,counter);
	
	// Create bounding sphere using Ritter [https://en.wikipedia.org/wiki/Bounding_sphere#cite_note-Ritter1990-2]
	//  Ritter, Jack (1990), "An efficient bounding sphere", in Glassner, Andrew S., Graphics Gems, San Diego, CA, US: Academic Press Professional, Inc., pp. 301–303, ISBN 0-12-286166-3
	double max_dist=0;
	
	Coords B_sphere_x = v1[0];
	Coords B_sphere_y = B_sphere_x ;
	for (i = 1 ; i < counter ; i++ ){
		if (sqrt(v1[0].x-v1[i].x+v1[0].y-v1[i].y+v1[0].z-v1[i].z) > max_dist){
			max_dist = sqrt(v1[0].x-v1[i].x+v1[0].y-v1[i].y+v1[0].z-v1[i].z);
			B_sphere_y = v1[i];
		}
		if (sqrt(v1[0].x-v2[i].x+v1[0].y-v2[i].y+v1[0].z-v2[i].z) > max_dist){
			max_dist = sqrt(v1[0].x-v2[i].x+v1[0].y-v2[i].y+v1[0].z-v2[i].z);
			B_sphere_y = v2[i];
		}
		if (sqrt(v1[0].x-v3[i].x+v1[0].y-v3[i].y+v1[0].z-v3[i].z) > max_dist){
			max_dist = sqrt(v1[0].x-v3[i].x+v1[0].y-v3[i].y+v1[0].z-v3[i].z);
			B_sphere_y = v3[i];
		}
	}
	Coords B_sphere_z = B_sphere_y ;
	max_dist=0;
	for (i = 1 ; i < counter ; i++ ){
		if (sqrt(pow(B_sphere_y.x-v1[i].x,2) + pow(B_sphere_y.y-v1[i].y,2) + pow(B_sphere_y.z-v1[i].z,2)) > max_dist){
			max_dist = sqrt(pow(B_sphere_y.x-v1[i].x,2) + pow(B_sphere_y.y-v1[i].y,2) + pow(B_sphere_y.z-v1[i].z,2));
			B_sphere_z = v1[i];
		}
		if (sqrt(pow(B_sphere_y.x-v2[i].x,2) + pow(B_sphere_y.y-v2[i].y,2) + pow(B_sphere_y.z-v2[i].z,2)) > max_dist){
			max_dist = sqrt(pow(B_sphere_y.x-v2[i].x,2) + pow(B_sphere_y.y-v2[i].y,2) + pow(B_sphere_y.z-v2[i].z,2));
			B_sphere_z = v2[i];
		}
		if (sqrt(pow(B_sphere_y.x-v3[i].x,2) + pow(B_sphere_y.y-v3[i].y,2) + pow(B_sphere_y.z-v3[i].z,2)) > max_dist){
			max_dist = sqrt(pow(B_sphere_y.x-v3[i].x,2) + pow(B_sphere_y.y-v3[i].y,2) + pow(B_sphere_y.z-v3[i].z,2));
			B_sphere_z = v3[i];
		}
	}
	
	double tmp_Radius = sqrt(pow(B_sphere_y.x-B_sphere_z.x,2.0) + pow(B_sphere_y.y-B_sphere_z.y,2.0) + pow(B_sphere_y.z-B_sphere_z.z,2.0))/2;
	Coords bbcenter = coords_set((B_sphere_y.x+B_sphere_z.x)/2 , (B_sphere_y.y+B_sphere_z.y)/2 , (B_sphere_y.z+B_sphere_z.z)/2);
	this_mesh_storage.Bounding_Box_Center = bbcenter;
	for (i = 0 ; i < counter ; i++ ){
		if (sqrt(pow(bbcenter.x-v1[i].x,2.0) + pow(bbcenter.y-v1[i].y,2.0) + pow(bbcenter.z-v1[i].z,2.0)) > tmp_Radius){
			tmp_Radius = sqrt(pow(bbcenter.x-v1[i].x,2.0) + pow(bbcenter.y-v1[i].y,2.0) + pow(bbcenter.z-v1[i].z,2.0));
		}
		if (sqrt(pow(bbcenter.x-v2[i].x,2.0) + pow(bbcenter.y-v2[i].y,2.0) + pow(bbcenter.z-v2[i].z,2.0)) > tmp_Radius){
			tmp_Radius = sqrt(pow(bbcenter.x-v2[i].x,2.0) + pow(bbcenter.y-v2[i].y,2.0) + pow(bbcenter.z-v2[i].z,2.0));
		}
		if (sqrt(pow(bbcenter.x-v3[i].x,2.0) + pow(bbcenter.y-v3[i].y,2.0) + pow(bbcenter.z-v3[i].z,2.0)) > tmp_Radius){
			tmp_Radius = sqrt(pow(bbcenter.x-v3[i].x,2.0) + pow(bbcenter.y-v3[i].y,2.0) + pow(bbcenter.z-v3[i].z,2.0));
		}
	
	}
	this_mesh_storage.Bounding_Box_Radius = tmp_Radius;
	printf("\n__________________________________\n BOUNDING BOX: \n\t Center = [%f , %f , %f]\n\t Radius = %f\n__________________________________\n", this_mesh_storage.Bounding_Box_Center.x,this_mesh_storage.Bounding_Box_Center.y,this_mesh_storage.Bounding_Box_Center.z,this_mesh_storage.Bounding_Box_Radius);
	
	
	
	
	// Convert to double arrays for further parsing (Coord has problems here - might be possible to do much faster!!)
	// Seriously not optimal... plz fix
	double v1_x[counter];
	double v1_y[counter];
	double v1_z[counter];
	double v2_x[counter];
	double v2_y[counter];
	double v2_z[counter];
	double v3_x[counter];
	double v3_y[counter];
	double v3_z[counter];
	double normal_x[counter];
	double normal_y[counter];
	double normal_z[counter];
	


	for (i = 0 ; i < counter ; i++ ){
		/*
		v1_x[i] = v1[i].x;
		v1_y[i] = v1[i].y;
		v1_z[i] = v1[i].z;
		
		v2_x[i] = v2[i].x;
		v2_y[i] = v2[i].y;
		v2_z[i] = v2[i].z;
		
		v3_x[i] = v3[i].x;
		v3_y[i] = v3[i].y;
		v3_z[i] = v3[i].z;
		
		normal_x[i] = normal[i].x;
		normal_y[i] = normal[i].y;
		normal_z[i] = normal[i].z;
		*/
		this_mesh_storage.v1_x[i] = v1[i].x;
		this_mesh_storage.v1_y[i] = v1[i].y;
		this_mesh_storage.v1_z[i] = v1[i].z;
		
		this_mesh_storage.v2_x[i] = v2[i].x;
		this_mesh_storage.v2_y[i] = v2[i].y;
		this_mesh_storage.v2_z[i] = v2[i].z;
		
		this_mesh_storage.v3_x[i] = v3[i].x;
		this_mesh_storage.v3_y[i] = v3[i].y;
		this_mesh_storage.v3_z[i] = v3[i].z;
		
		this_mesh_storage.normal_x[i] = normal[i].x;
		this_mesh_storage.normal_y[i] = normal[i].y;
		this_mesh_storage.normal_z[i] = normal[i].z;
		
		//printf("\n v2 [%i]: v2[%i] = [%f,%f,%f] \n",iter,v2_x[i],v2_y[i],v2_z[i]);
	}

// check if this is a cylinder and add coordinates of tip if it is a mesh

this_mesh_storage.counter = counter;
this_mesh_storage.n_facets = n_facets;


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

this_mesh_volume.geometry.geometry_p_interact = p_interact;



this_mesh_volume.geometry.visualization_on = visualize;
this_mesh_volume.geometry.geometry_parameters.p_mesh_storage = &this_mesh_storage;
this_mesh_volume.geometry.within_function = &r_within_mesh;
this_mesh_volume.geometry.intersect_function = &sample_mesh_intersect;
this_mesh_volume.geometry.mcdisplay_function = &mcdisplay_mesh_function;
this_mesh_volume.geometry.shell_points = &mesh_shell_points;
this_mesh_volume.geometry.initialize_from_main_function = &initialize_mesh_geometry_from_main_component;
this_mesh_volume.geometry.process_rot_allocated = 0;
this_mesh_volume.geometry.copy_geometry_parameters = &allocate_mesh_storage_copy;
rot_copy(this_mesh_volume.geometry.rotation_matrix,ROT_A_CURRENT_COMP);
rot_transpose(ROT_A_CURRENT_COMP,this_mesh_volume.geometry.transpose_rotation_matrix);


// Initialize loggers
this_mesh_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,"%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_mesh_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
%{
%}

END