/*******************************************************************************
*         McXtrace instrument definition URL=http://www.mcxtrace.org
*
* Instrument: SOLEIL_CRISTAL
*
* %Identification
* Written by: FARHI Emmanuel with help of Claude (Anthropic), adapted from SOLEIL_MARS and SOLEIL_DIFFABS
* Date: September 2026
* Origin: SOLEIL
* Version: 0.1
* %INSTRUMENT_SITE: SOLEIL
*
* SOLEIL CRISTAL beamline: diffraction with 2/4/6-circle diffractometers.
*
* %Description
* CRISTAL is a diffraction beam-line covering 5-30 keV (CRISTAL_dE/E ~1e-4), fed by an
* in-vacuum U20 undulator (98 periods, 20 mm period). The optical layout, in
* beam order, is: primary slit, Si(111) double-crystal monochromator (DCM,
* second crystal sagittally bent for horizontal focusing - not modelled here
* for simplicity), a pair of horizontally-deflecting Si mirrors (Rh/Pt coated,
* first one bendable, for horizontal focusing and harmonic rejection), a
* secondary slit, and finally the sample stage. Most powder-diffraction work
* is done on the 2-circle diffractometer, equipped with a curved Mythen-II
* pixel detector (radius 720 mm). Only components that act on the beam are
* modelled; hutch diagnostics such as xbpm, beam imagers and attenuators are
* omitted for simplicity.
*
* %Example: CRISTAL_E0=12 Detector: CRISTAL_detector_diffraction_I=700000.0
*
* %Parameters
* CRISTAL_E0:        [keV] Central energy of the beam
* CRISTAL_dE:        [keV] Energy half-bandwidth at the photon source
* CRISTAL_M1_angle:  [mrad] Grazing angle of the first (bendable) horizontal mirror
* CRISTAL_M2_angle:  [mrad] Grazing angle of the second horizontal mirror
* CRISTAL_M1_radius: [m]   Bending radius of mirror M1 (0=flat)
* CRISTAL_M2_radius: [m]   Bending radius of mirror M2 (0=flat)
* CRISTAL_M_coating: [str] Mirror coating (Rh or Pt stripe), e.g. "Rh.txt" or "Pt.dat"
* reflections:[str] Sample structure file, LAU/CIF format
*
* %Link
* https://www.synchrotron-soleil.fr/en/beamlines/cristal
*
* %End
*******************************************************************************/

DEFINE INSTRUMENT SOLEIL_CRISTAL(CRISTAL_E0=12, CRISTAL_dE=1e-4,
  CRISTAL_M1_angle=3, CRISTAL_M2_angle=3, CRISTAL_M1_radius=0, CRISTAL_M2_radius=0,
  string CRISTAL_M_coating="Rh.txt",
  string CRISTAL_reflections="LaB6_660b_AVID2.hkl",
  int splits=10)

DECLARE
%{
  // Primary slit (Monochromator Ring Hutch, before the DCM)
  double CRISTAL_DistanceSourceToSlit1 = 14.0;
  double CRISTAL_Slit1Width  = 0.02;
  double CRISTAL_Slit1Height = 0.01;

  // Si(111) DCM
  double CRISTAL_DistanceSourceToDCM = 17.0;
  double CRISTAL_dcm_gap = 0.01;
  double CRISTAL_DCM_theta;

  // Mirror pair (Mirrors Hutch, horizontal deflection)
  double CRISTAL_DistanceSourceToM1 = 20.0;
  double CRISTAL_DistanceM1ToM2     = 0.5;
  double CRISTAL_M_length = 0.3;
  double CRISTAL_M_width  = 0.1;

  // Secondary slit (end of Mirrors Hutch)
  double CRISTAL_DistanceSourceToSlit2 = 22.5;
  double CRISTAL_Slit2Width  = 0.03;
  double CRISTAL_Slit2Height = 0.02;

  // Sample stage (2-circle powder diffractometer, BC Hutch)
  double CRISTAL_DistanceSourceToSample = 35.0;

  // Mythen-II curved detector, radius 720 mm
  double CRISTAL_DetectorRadius = 0.72;
%}

INITIALIZE
%{
  // Si(111): a=5.4309 Angs, d111 = a/sqrt(3)
  {
    double d = 5.4309/sqrt(3);
    double sin_theta = 12.398/(2*d*CRISTAL_E0);

    if (fabs(sin_theta) >= 1 || CRISTAL_E0 < 5 || CRISTAL_E0 > 30)
      exit(fprintf(stderr,
        "%s: ERROR: Si(111) can not reflect CRISTAL_E0=%g keV (valid range ~5-30 keV). Aborting.\n",
        NAME_INSTRUMENT, CRISTAL_E0));

    CRISTAL_DCM_theta = asin(sin_theta)*RAD2DEG;
  }

  MPI_MASTER(
  printf("%s: CRISTAL_E0=%g [keV] CRISTAL_dE=%g [keV]\n", NAME_INSTRUMENT, CRISTAL_E0, CRISTAL_dE);
  printf("%s: DCM theta=%g [deg]\n",   NAME_INSTRUMENT, CRISTAL_DCM_theta);
  printf("%s: Mirrors M1=%g [mrad] M2=%g [mrad]\n", NAME_INSTRUMENT, CRISTAL_M1_angle, CRISTAL_M2_angle);
  );

  CRISTAL_M1_angle *= RAD2DEG/1000;
  CRISTAL_M2_angle *= RAD2DEG/1000;
%}

TRACE

// electronic beam spec along storage ring 354m diameter
// BeamLineName| Long_Pos(m)| RMS_H_Size(µm)| RMS_V_Size(µm)| RMS_H_div(µrad)| RMS_V_div(µrad)
// -------|--------|-----|----|----|--
// CRISTAL|121,7489|332,9|9,3|17,3|3,6
REMOVABLE COMPONENT CRISTAL_origin = Progress_bar()
AT (0, 0, 0) RELATIVE ABSOLUTE

/* -------------------------------------- Source: in-vacuum U20 undulator */
// Undulator (in-vacuum, 98 periods of 20 mm)
COMPONENT CRISTAL_Source = Undulator(
  E0 = CRISTAL_E0,
  dE = CRISTAL_dE,
  Ee = 2.75,
  Ie = 0.5,
  K = 5,
  sigex = 388e-6,
  sigey = 8.1e-6,
  sigepx = 14.5e-6,
  sigepy = 4.61e-6,
  focus_xw=CRISTAL_Slit1Width, focus_yh=CRISTAL_Slit1Height,
  dist=CRISTAL_DistanceSourceToSlit1)
AT (0, 0, 0) RELATIVE PREVIOUS

/* -------------------------------------- Primary slit s1 */
COMPONENT CRISTAL_slit1 = Slit(
    xwidth=CRISTAL_Slit1Width, yheight=CRISTAL_Slit1Height)
AT (0, 0, CRISTAL_DistanceSourceToSlit1) RELATIVE CRISTAL_Source

/* -------------------------------------- Si(111) DCM */
COMPONENT CRISTAL_dcm_location = Arm()
AT (0, 0, CRISTAL_DistanceSourceToDCM-CRISTAL_DistanceSourceToSlit1) RELATIVE PREVIOUS

COMPONENT CRISTAL_dcm_xtal0 = Bragg_crystal(
    length=0.1, width=0.05,
    alpha=0, h=1, k=1, l=1, material="Si.txt", crystal_type=2)
AT (0, 0, 0)             RELATIVE CRISTAL_dcm_location
ROTATED (-CRISTAL_DCM_theta,0,0) RELATIVE CRISTAL_dcm_location
EXTEND
%{
  if (!SCATTERED) ABSORB;
%}

COMPONENT CRISTAL_arm_xtal0 = Arm()
AT (0, 0, 0)             RELATIVE PREVIOUS
ROTATED (-CRISTAL_DCM_theta,0,0) RELATIVE PREVIOUS

// second crystal: kept flat here for simplicity (sagittal bending for
// horizontal focusing, as used on CRISTAL, is not modelled)
COMPONENT CRISTAL_dcm_xtal1 = COPY(CRISTAL_dcm_xtal0)
AT (0, CRISTAL_dcm_gap, CRISTAL_DCM_theta ? CRISTAL_dcm_gap/tan(CRISTAL_DCM_theta*DEG2RAD) : 0) RELATIVE CRISTAL_dcm_xtal0
ROTATED (CRISTAL_DCM_theta,0,0) RELATIVE CRISTAL_arm_xtal0
EXTEND
%{
  if (!SCATTERED) ABSORB;
%}

COMPONENT CRISTAL_arm_xtal1 = Arm()
AT (0, 0, 0)            RELATIVE PREVIOUS
ROTATED (CRISTAL_DCM_theta,0,0) RELATIVE PREVIOUS

/*Diagnostic: energy after the DCM*/
COMPONENT CRISTAL_emon_dcm = E_monitor(
    xwidth=0.05, yheight=0.05, filename="CRISTAL_emon_dcm",
    Emin=CRISTAL_E0-3*CRISTAL_dE, Emax=CRISTAL_E0+3*CRISTAL_dE, nE=101)
AT (0, 0, 0.05) RELATIVE CRISTAL_arm_xtal1

/* -------------------------------------- Mirror M1 (bendable, horizontal) */
COMPONENT CRISTAL_M1_location = Arm()
AT (0, 0, CRISTAL_DistanceSourceToM1-CRISTAL_DistanceSourceToDCM-0.05) RELATIVE PREVIOUS

COMPONENT CRISTAL_M1_rotated = Arm()
AT (0, 0, 0)          RELATIVE CRISTAL_M1_location
ROTATED (-CRISTAL_M1_angle,0,0) RELATIVE CRISTAL_M1_location

COMPONENT CRISTAL_M1 = Mirror_curved(
    length=CRISTAL_M_length, width=CRISTAL_M_width,
    coating=CRISTAL_M_coating, radius=CRISTAL_M1_radius)
AT (0, 0, 0)       RELATIVE CRISTAL_M1_rotated
ROTATED (0, 0, 90) RELATIVE CRISTAL_M1_rotated
EXTEND
%{
  if (!SCATTERED) ABSORB;
%}

COMPONENT CRISTAL_M1_out = Arm()
AT (0, 0, 0)             RELATIVE CRISTAL_M1_rotated
ROTATED (-CRISTAL_M1_angle,0,0)  RELATIVE CRISTAL_M1_rotated

/* -------------------------------------- Mirror M2 (fixed, horizontal) */
COMPONENT CRISTAL_M2_location = Arm()
AT (0, 0, CRISTAL_DistanceM1ToM2) RELATIVE CRISTAL_M1_out

COMPONENT CRISTAL_M2_rotated = Arm()
AT (0, 0, 0)           RELATIVE CRISTAL_M2_location
ROTATED (CRISTAL_M2_angle,0,0) RELATIVE CRISTAL_M2_location

COMPONENT CRISTAL_M2 = Mirror_curved(
    length=CRISTAL_M_length, width=CRISTAL_M_width,
    coating=CRISTAL_M_coating, radius=CRISTAL_M2_radius)
AT (0, 0, 0)        RELATIVE CRISTAL_M2_rotated
ROTATED (0, 0, -90) RELATIVE CRISTAL_M2_rotated
EXTEND
%{
  if (!SCATTERED) ABSORB;
%}

COMPONENT CRISTAL_M2_out = Arm()
AT (0, 0, 0)           RELATIVE CRISTAL_M2_rotated
ROTATED (CRISTAL_M2_angle,0,0) RELATIVE CRISTAL_M2_rotated

/* -------------------------------------- Secondary slit s2 */
COMPONENT CRISTAL_slit2 = Slit(
    xwidth=CRISTAL_Slit2Width, yheight=CRISTAL_Slit2Height)
AT (0, 0, CRISTAL_DistanceSourceToSlit2-CRISTAL_DistanceSourceToM1-CRISTAL_DistanceM1ToM2) RELATIVE PREVIOUS

/*Diagnostic: beam profile after the mirrors*/
COMPONENT CRISTAL_psd_slit2 = PSD_monitor(filename="CRISTAL_psd_slit2")
AT (0, 0, 0) RELATIVE PREVIOUS

/* -------------------------------------- Sample stage (2-circle diffractometer) */
COMPONENT CRISTAL_sample_stage = Arm()
AT (0, 0, CRISTAL_DistanceSourceToSample-CRISTAL_DistanceSourceToSlit2) RELATIVE PREVIOUS

// The Rayleigh scattering at CRISTAL_E0 is large and creates massive background;
// it is ignored here, as in the MARS/DIFFABS models.
SPLIT splits COMPONENT CRISTAL_fluo = FluoPowder(
    material=CRISTAL_reflections,
    xwidth=0.001, yheight=0.01, zdepth=5e-6,
    focus_aw=180, focus_ah=asin(0.02/CRISTAL_DetectorRadius)*RAD2DEG, target_z=1, flag_rayleigh=0,
    p_interact=0.99, d_phi=asin(0.02/CRISTAL_DetectorRadius)*RAD2DEG)
AT (0, 0, 0) RELATIVE CRISTAL_sample_stage

/* -------------------------------------- Detector: curved Mythen-II, R=720 mm */
COMPONENT CRISTAL_detector_diffraction = Monitor_nD(
    bins=50000, options="abs theta", min=2, max=100,
    radius=CRISTAL_DetectorRadius, yheight=0.02, restore_xray=1)
WHEN (0.75*CRISTAL_E0 < K2E*sqrt(kx*kx+ky*ky+kz*kz)) // model detector energy discrimination
AT (0, 0, 0) RELATIVE CRISTAL_sample_stage

COMPONENT CRISTAL_detector_fluo = Monitor_nD(
    bins=1024, options="energy", min=0, max=CRISTAL_E0*1.2,
    radius=CRISTAL_DetectorRadius, yheight=0.02)
WHEN (K2E*sqrt(kx*kx+ky*ky+kz*kz) < 0.95*CRISTAL_E0)
AT (0, 0, 0) RELATIVE CRISTAL_sample_stage

END