/*******************************************************************************
* Instrument: ROCK beam-line at SOLEIL
*
* %Identification
* Written by: Stephane Bac, Antoine Padovani, Emmanuel Farhi
* Date: 23/02/2022
* Origin: SOLEIL
* Version: 0.3
* %INSTRUMENT_SITE: SOLEIL
*
* ROCK beam-line at SOLEIL
*
* %Description
* ROCK : Rocking Optics for Chemical Kinetics ROCK time-resolved X-ray
*
* Absorption spectroscopy (XAS) beamline Energy range 4 - 40 keV 
* The ROCK beamline (ROCK being the acronym for Rocking Optics for Chemical Kinetics)
* is devoted to the study of fast kinetic processes in nanomaterials used in
* catalysis and batteries. The objective is to contribute to the development of
* more efficient catalysts and batteries which should find successful industrial
* applications in the field of energy generation and storage in compliance with
* the protection of public health and environment. The better knowledge at the
* atomic scale of nanomaterials involved in catalysis or energy storage provided
* by time-resolved XAS is recognized by the concerned communities as mandatory
* for establishing synthesis strategies leading to important breakthroughs in
* the production of energy from renewable sources and in the development of
* advanced energy storage devices. 
*
* Position | Element
* ---------|--------------------------------------------------------------------
* 0        | Bending_magnet 1.72 T
* 8.5336   | Slit 1/2
* 10.15    | Toroidal mirror M1
* 11.69    | Slit 3
* 16.82    | Mirror M2a
* 18.15    | Slit 4
* 19       | Channel-cut Si monochromator at 20 (Si220), 18.92 or 19.25 m (Si111)
* 21.13    | Slit 5/6
* 22.44    | Mirror M2b
* 32.44    | Sample
*
* Examples:
* Copper scan: mxrun SOLEIL_ROCK.instr E0=8.500,9.500 scan=1 cc=2 -N100
* Manganese and chrome scan: mxrun SOLEIL_ROCK.instr E0=5.700,6.800 scan=1 sample_file=MnCr cc=2 -N100
* Pretty energy repartition monitor result: mxrun SOLEIL_ROCK.instr E0=17.000 cc=2 -n1e8
*
* %Example: E0=15.918 Detector: psd_monitor_I=6.9e+08
*
* %Parameters
* E0:               [keV]    Energy to hit, i.e. selected by the channel-cut monochromator
* dE:               [keV]    Energy spread at the source
* cc:               [0-3]    Channel-cut monochromator type. 0 for Si 220 with an energy range of  5.62883-46.2834 eV/4-35 deg (hit-table:11.752-34.055 ev/5.44-15.94 deg); 1 for Si 111 long 3.44694-28.3427 eV/4-35 deg (hit-table:7.196-20.854 ev/5.44-15.94 deg); 2 for Si 111 short 3.44694-18.914.3 eV/6-35 deg (hittable:5.323-18.914 ev/6-21.8 deg); 3 changes cc dynamically
* scan:             [0-1]    0 no energy scan, 1 energy scan
* angle_m2a_m2b:    [rad]    M2A/M2B mirror's deviation angle, can vary from 0.0035 to 0.0104
* angle_m1:         [rad]    M1 mirror's deviation angle
* sample_file:      [string] Sample chemical formulae
* reflec_material_M1:      [str] Material reflectivity file name for M1 mirror, e.g. "Ir.dat"
* reflec_material_M2A_M2B: [str] Material reflectivity file name for M2A and M2B mirror. Use NULL for automatic setting.
* 
* %Link
* https://www.synchrotron-soleil.fr/en/beamlines/rock
* %Link
* https://gitlab.synchrotron-soleil.fr/grades/mcxtrace-rock
* %Link
* https://github.com/antoinepado/glitch_runner_rock
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT SOLEIL_ROCK (E0=15.918, dE=1,
    scan=0, 
    angle_m2a_m2b=0.008, 
    angle_m1=0.0045,
    string sample_file="CuMo", 
    string reflec_material_M1="Ir.dat",
    string reflec_material_M2A_M2B="NULL",
    int cc=0)

DECLARE
%{
    double calculated_angle;
    double Eminkev;
    double Emaxkev;
    double incr;
    int n;
    int i_h;
    int i_k;
    int i_l;
    double length_first_crystal;
    double length_second_crystal;
    double center_first_crystal;
    double dist_center_to_center;
    double center_to_c_mono_angle;
    char e_repartition_options[512];     
    double I_err_calc; 
    double absorption_coefficient;
    double absorption_coefficient_error;   
    double number_events_after_sample;       
%}

INITIALIZE
%{
    i_h=1;
    i_k=1;
    i_l=1;
    
    if(cc>=3){ //automatically select the right cc for the energy E0
        if(E0<=15.516){ 
            cc = 2;
        }
        else if (E0<=17.439){
            cc = 1;
        }
        else if (E0<=34.055){
            cc = 0;	
        }
    }

    if(cc<=0){ //cc 220
        Eminkev=5.62883;
        Emaxkev=46.2834;
        incr=0.995/(Emaxkev-Eminkev);
        n=8;
        i_h=2;
        i_k=2;
        i_l=0;
        length_first_crystal=70e-3;
        length_second_crystal=70e-3;
        center_first_crystal=20;
        dist_center_to_center = sqrt(0.01*0.01+0.07*0.07);
        /// center_to_c_mono_angle is arctan(e/(length_first_crystal/2 + length_second_crystal/2))
        center_to_c_mono_angle=0.14189;
    }
    else if (cc<=1){ //cc long 111
        Eminkev=3.44694;
        Emaxkev=28.3427;
        incr=0.995/(Emaxkev-Eminkev);
        n=3;
        length_first_crystal=70e-3;
        length_second_crystal=70e-3;
        center_first_crystal=18.92;
        dist_center_to_center = sqrt(0.01*0.01+0.07*0.07);
        center_to_c_mono_angle=0.14189;
    } else if (cc<=2){ //cc short 111
        Eminkev=3.44694;
        Emaxkev=18.9143;
        incr=0.995/(Emaxkev-Eminkev);
        n=3;
        length_first_crystal=50e-3;
        length_second_crystal=70e-3;
        center_first_crystal=19.25;
        dist_center_to_center = sqrt(0.01*0.01+0.06*0.06);
        center_to_c_mono_angle=0.1651;
    }

    calculated_angle=RAD2DEG*asin(12.39842*sqrt(n)/(2*5.4309*E0));
    
    if (!reflec_material_M2A_M2B || !strlen(reflec_material_M2A_M2B) 
        || !strcmp(reflec_material_M2A_M2B, "NULL") || !strcmp(reflec_material_M2A_M2B, "0")) {

      if(E0<=7.98){
          reflec_material_M2A_M2B="B4C.dat";
      }
      else if (E0<=16.46){
          reflec_material_M2A_M2B="Pd.dat";
      }
      else if (E0>16.46){
          reflec_material_M2A_M2B="Pt.dat";
      }
    }
    
    MPI_MASTER(
    fprintf(stdout,"%s: Energy=%g [keV]\n"
      "\tM2 mirror=%s\n"
      "\tMonochromator Si<%i %i %i> angle %g [deg], config=%i\n", 
      NAME_INSTRUMENT, E0, reflec_material_M2A_M2B, i_h,i_k,i_l, calculated_angle, cc);
    );

    sprintf(e_repartition_options,"energy limits=[%g %g], y",
        E0-0.02*E0, E0+0.02*E0);  
%}

TRACE

COMPONENT Origin = Progress_bar(
) 
AT (0, 0, 0) ABSOLUTE

/* -------------------------------------------------- Source */
COMPONENT Source = Bending_magnet(
    E0=E0, dE = dE,
    Ee = 2.75, Ie = 0.5, B = 1.72, sigex=54.9e-6, sigey=20.2e-6
    )
AT (0, 0, 0) RELATIVE Origin

/* -------------------------------------------------- Vertical slit 1 */
COMPONENT slit1 = Slit(
    xwidth=0.021, 
    yheight=0.042)
AT (0, 0, 8.5336) RELATIVE PREVIOUS

/* -------------------------------------------------- Horizontal slit 2 */
COMPONENT slit2 = Slit(
    xwidth=0.02, 
    yheight=0.01)
AT (0, 0, 8.6486-8.5336) RELATIVE PREVIOUS

/* -------------------------------------------------- Toroidal mirror M1 */   
COMPONENT mirror_m1 = Mirror_toroid_pothole(
    xwidth = 0.015,
    zdepth = 1.1,
    radius=0.0317, 
    radius_o=9.02e3,coating=reflec_material_M1)    
AT (0, 0, 10.15-8.6486) RELATIVE PREVIOUS	
ROTATED (-angle_m1*RAD2DEG/2, 0, 0) RELATIVE PREVIOUS
EXTEND
%{ 
	if (!SCATTERED) ABSORB;
%}

COMPONENT mirror_m1_out = Arm()
AT (0,0,0) RELATIVE PREVIOUS
ROTATED (-angle_m1*RAD2DEG/2, 0, 0) RELATIVE PREVIOUS

/* -------------------------------------------------- Horizontal slit  3 */

COMPONENT slit3_xy = PSD_monitor(xwidth=0.02, yheight=0.01)
AT (0, 0, 11.6905-10.15) RELATIVE PREVIOUS

COMPONENT slit3_div = Divergence_monitor(maxdiv_h=0.05, maxdiv_v=0.05, nh=128, nv=128)
AT (0, 0, 0) RELATIVE PREVIOUS

COMPONENT slit3 = Slit(
    xwidth=0.02, 
    yheight=0.01)
AT (0, 0, 0) RELATIVE slit3_xy

/* -------------------------------------------------- M2a mirror */
COMPONENT mirror_m2a = Mirror(
    xwidth=47e-3,
    zdepth=1100e-3,coating=reflec_material_M2A_M2B)
AT (0, 0, (16.82-11.6905)) RELATIVE PREVIOUS	
ROTATED (-angle_m2a_m2b*RAD2DEG/2, 0, 0) RELATIVE PREVIOUS
EXTEND
%{ 
	if (!SCATTERED) ABSORB;
%}

COMPONENT mirror_m2a_out = Arm()
AT (0, 0, 0) RELATIVE PREVIOUS
ROTATED (-angle_m2a_m2b*RAD2DEG/2, 0, 0) RELATIVE PREVIOUS

/* -------------------------------------------------- Horizontal slit 4 */

COMPONENT slit4_in = PSD_monitor(xwidth=0.5, yheight=0.5)
AT (0, 0, 18.1512-16.82) RELATIVE mirror_m2a_out

COMPONENT slit4 = Slit(
    xwidth=0.02, 
    yheight=0.02)
AT (0, 0, 0) RELATIVE PREVIOUS	

/* -------------------------------------------------- CC crystal 1 */

COMPONENT arm_crystal0 = Arm()
AT (0, 0, center_first_crystal-18.1512) RELATIVE slit4

COMPONENT bragg_crystal = Bragg_crystal(
    length=length_first_crystal, 
    width=25e-3,
    material = "Si.txt",
    h = i_h,
    k = i_k,
    l = i_l,
    crystal_type = 2
)
AT (0, 0, 0) RELATIVE arm_crystal0	
ROTATED (-calculated_angle, 0, 0) RELATIVE arm_crystal0
EXTEND
%{ 
	if (!SCATTERED) ABSORB;
%}

COMPONENT arm_crystal1 = Arm()
AT (0, 0, 0) RELATIVE PREVIOUS
ROTATED (-calculated_angle, 0, 0) RELATIVE PREVIOUS

/* -------------------------------------------------- CC crystal 2 */
// gap btw crystals is 1 cm. We shit the 2nd blade so that it is centred wrt reflection.
COMPONENT bragg_crystal_two = COPY(bragg_crystal)(
    length=length_second_crystal) 
AT (0, 0.01, calculated_angle ? 0.01/tan(calculated_angle*DEG2RAD) : 0) RELATIVE bragg_crystal
ROTATED (calculated_angle, 0, 0) RELATIVE arm_crystal1
EXTEND
%{ 
	if (!SCATTERED) ABSORB;
%}

COMPONENT arm_crystal2 = Arm()
AT (0, 0, 0) RELATIVE PREVIOUS
ROTATED (calculated_angle, 0, 0) RELATIVE PREVIOUS

COMPONENT bragg_crystal_out = PSD_monitor(xwidth=0.1, yheight=0.1, nx=200, ny=200)
AT (0,0,2*length_second_crystal)  RELATIVE arm_crystal2

/* -------------------------------------------------- Horizontal slit 5 */
// the beam is positioned at Y=+dy*sin(2*calculated_angle*DEG2RAD)/sin(calculated_angle*DEG2RAD) 
// with dy=0.01 (CC gap) wrt the CC entry
COMPONENT slit5 = Slit(
    xwidth=0.02, 
    yheight=0.01)
AT (0, 0, 21.135 - center_first_crystal) RELATIVE arm_crystal2


/* -------------------------------------------------- Vertical slit 6 */
COMPONENT slit6 = Slit(
    xwidth=0.01, 
    yheight=0.02)
AT (0, 0, 21.25-21.135) RELATIVE PREVIOUS

/* -------------------------------------------------- M2b mirror */    
COMPONENT mirror_m2b = Mirror_curved(
    radius=5e3, 
    width=47e-3,
    length=1100e-3,
    coating=reflec_material_M2A_M2B)          
AT (0, 0, 22.44-21.25) RELATIVE PREVIOUS	
ROTATED (angle_m2a_m2b*RAD2DEG/2, 0, 90) RELATIVE PREVIOUS
EXTEND
%{ 
	if (!SCATTERED) ABSORB; 
%}

COMPONENT mirror_m2b_out = Arm()
AT (0, 0, 0) RELATIVE PREVIOUS
ROTATED (angle_m2a_m2b*RAD2DEG, 0, 0) RELATIVE arm_crystal2 

/* -------------------------------------------------- E monitor */ 
COMPONENT EnergyMonitor_before_sample = Monitor_nD(  
    xwidth =0.1 , 
    yheight =0.1 ,
    min=E0-(0.005+(E0-Eminkev)*incr), 
    max=E0+(0.005+(E0-Eminkev)*incr), 
    bins=500,
    options="energy",
    filename="EnergyMonitor_before_sample",
    restore_xray = 1
    ) 
AT (0,0,10) RELATIVE PREVIOUS



/* -------------------------------------------------- sample */
// this is the sample stage
COMPONENT e_repartition_after_m2b = Monitor_nD(
  xwidth =0.02 , yheight =0.02, 
  options=e_repartition_options, bins=500,
  filename="e_repartition_after_m2b",
  restore_xray = 1) 
AT (0,0,0.5) RELATIVE PREVIOUS

COMPONENT sample_in = PSD_monitor(xwidth=0.01, yheight=0.01, nx=200, ny=200)
AT (0,0,0) RELATIVE PREVIOUS

SPLIT COMPONENT absorption_sample = Fluorescence(
  material=sample_file,
  xwidth = 0.05,
  yheight = 0.05,
  zdepth = 0.0001)
AT (0, 0, 0) RELATIVE e_repartition_after_m2b

COMPONENT fluo_monitor = Monitor_nD(
  radius=0.2, options="energy", bins=1024, min=0, max=E0*1.2)
AT (0,0,0) RELATIVE PREVIOUS

/* -------------------------------------------------- E monitor */
COMPONENT EnergyMonitor_after_sample = COPY(EnergyMonitor_before_sample)(
        filename     = "EnergyMonitor_after_sample"
)  
AT (0,0,0.5) RELATIVE PREVIOUS

COMPONENT psd_monitor = PSD_monitor(
    filename="psd", 
    xwidth=0.004, 
    yheight=0.004,            
    nx=100,
    ny=100,
    restore_xray = 1
    )
AT (0, 0, 0) RELATIVE PREVIOUS


FINALLY
%{
    // If there is a scan, calculate the absorption coefficient.
    if(scan>=1){ 
#ifdef USE_MPI
        //MPI_MASTER is equivalent to if(mpi_node_root>=mpi_node_rank){ //statement }
        //If DETECTOR_OUT_0D is inside of MPI_MASTER the program hangs forever. 
        //If outside it seems to work. Weird.
        MPI_MASTER(
#endif 
        MCDETECTOR e_before_sample = COMP_GETPAR2(EnergyMonitor_before_sample,detector);
        MCDETECTOR e_after_sample =  COMP_GETPAR2(EnergyMonitor_after_sample,detector);
        absorption_coefficient = log(e_before_sample.intensity/e_after_sample.intensity);              
        I_err_calc = ((e_before_sample.error/e_before_sample.intensity)+(e_after_sample.error/e_after_sample.intensity));                  
        absorption_coefficient_error = I_err_calc;  
        number_events_after_sample = e_after_sample.events; 

        fprintf(stdout,"*** \n");    
        fprintf(stdout,"Before the sample: I %g I_error %g N %g \n", e_before_sample.intensity, e_before_sample.error, e_before_sample.events);
        fprintf(stdout,"After the sample: I %g I_error %g N %g \n", e_after_sample.intensity, e_after_sample.error, e_after_sample.events); 
        fprintf(stdout,"I_err_calc: I %g \n", I_err_calc);
        fprintf(stdout,"*** \n ");  
                                                 
#ifdef USE_MPI                  
        );
        //MPI_Barrier(MPI_COMM_WORLD);
#endif  
                      
        // This set of defines is to avoid getting a '.' in the component name
        #ifdef NAME_CURRENT_COMP
        #undef NAME_CURRENT_COMP
        #define NAME_CURRENT_COMP "absorption_coefficient"
        #endif           
	Coords dummy;
	Rotation Rot;
	rot_set_rotation(Rot,0,0,0);
	mcdetector_out_0D("XANES EXAFS",number_events_after_sample, absorption_coefficient, absorption_coefficient_error,instrument_name,dummy,Rot);  
    }
%}

END