/*******************************************************************************
*         McStas instrument definition URL=http://www.mcstas.org
*
* Instrument: Test_Incoherent_MS
*
* %Identification
* Written by: Daniel Lomholt Christensen and Peter Willendrup
* Date: 20250603
* Origin: NBI
* %INSTRUMENT_SITE: Tests_samples
*
* Test instrument for validation of magnitude of multiple scattering in Incoherent.comp.
*
* %Description
* Test instrument for validation of magnitude of multiple scattering in Incoherent.comp. 
* The instrument considers "thin" and "thick" sample conditions.
*
*
* For the thin sample conditions, the following analytical calulations find the 
* Neutron intensity at the detector to be 27.12:
* <div class="latex">
* $$ 
* \frac{dN_{\rm real}}{dt}= \Psi A \rho_{\rm c}\sigma_{\rm inc}
*                   \frac{\Delta \Omega}{4 \pi} a_{\rm lin} l_{unscattered}.
* $$
* Where $\Psi=1e7$, $A=1$, $\rho=0.0723$, $\sigma_{\rm inc}=\sigma_{\rm abs} = 5.08$, $\frac{\Delta \Omega}{4 \pi} = 10/99.95^2$ 
* and $a_{\rm lin}= \exp(-(\sigma_{\rm inc} + \sigma_{\rm abs}\lambda/1.7982)\rho l_{avg})$. 
* $l_{unscattered}=0.1cm$ and is the path through the sample eg. unscattered path length.
* </div>
*
*
* For the thick sample condition, the general formula is the same,
* But the attenuation is calculated using the effective length, instead of the actual
* length:
* <div class="latex">
* \begin{equation} \label{eq:inc_layer_l}
* l_{\rm app} = \int _{0}^{l_{\rm max}} a_{\rm lin}(z) dz 
* =  \int _{0}^{l_{\rm max}} e^{-2\mu z}dz 
* =  \frac{1-e^{-2\mu l_{\rm max}}}{2\mu },
* \end{equation}
* </div>
*
* By scanning across thickness values, one then finds a curve, that should fit the
* analytical calculations.
* 
*
* %Example: mcrun Test_Incoherent_MS sample=none Detector: total_scat_I=9.99909e+06
* %Example: mcrun Test_Incoherent_MS sample=thin Detector: psd_det_I=27.12
* %Example: mcrun Test_Incoherent_MS sample=thick Detector: psd_det_I=26.59
* %Example: mcrun Test_Incoherent_MS sample=thick Detector: psd_det_I=26.59
*
* %Parameters
* flux_mult: [1/(s*cm**2*st*meV)] Source flux x 1e14 multiplier for source 
* thick:                      [m] Sample thickness
* sample:                   [str] Sample state "thick", "thin", "none" supported
* total_scattering:           [1] Flag to indicate if sample focuses (=0) or not (=1)
* order:                      [1] Maximal order of multiple scattering
*
* %Link
* A reference/HTML link for more information
*
* %End
*******************************************************************************/
DEFINE INSTRUMENT Test_Incoherent_MS(flux_mult=0.5512,thick=0.001, 
				     int total_scattering=0, string sample="thin", int order=1)

/* The DECLARE section allows us to declare variables or  small      */
/* functions in C syntax. These may be used in the whole instrument. */
DECLARE
%{
  double det_rot;
  double E_i;
  double focus_w;
  double focus_h;
  int sample_case=0; // 0 = none, 1 = thin, 2 = thick
  #pragma acc declare create(sample_case)
%}

/* The INITIALIZE section is executed when the simulation starts     */
/* (C code). You may use them as component parameter values.         */
INITIALIZE
%{

  /* Set up default sample focusing */
  focus_w = 0.031822776602;
  focus_h = 0.031822776602;
  
  /* Determine if this is "thin", "thick" or "hollow" sample case */
  if (!strcmp("thin", sample)){
    det_rot=20;
    E_i = 25.3;
    sample_case=1;
  } else if (!strcmp("thick", sample)){
    det_rot=180;
    E_i = 10;
    sample_case=2;
  } else if (!strcmp("none", sample)){
    det_rot=20;
    E_i = 25.3;
    total_scattering=1;
  } else {
    fprintf(stderr,"ERROR: Sample case %s is not supported!\n",sample);
  }
  #pragma acc update device(sample_case)
  
  if (total_scattering){
    focus_w = 0;
    focus_h = 0;
  }

%}

/* Here comes the TRACE section, where the actual      */
/* instrument is defined as a sequence of components.  */
TRACE

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

COMPONENT source = Source_simple(
    yheight = 0.01,
    xwidth = 0.01,
    dist = 10.5,
    focus_xw = .01,
    focus_yh = .01,
    E0 = E_i,
    dE = 0.1,
    flux = 1e14*flux_mult
) AT (0, 0, 0) RELATIVE origin
ROTATED (0, 0, 0) RELATIVE origin

COMPONENT pre_samp_arm = Arm(
    ) AT (0, 0, 10.475) RELATIVE source

COMPONENT pre_samp_det = PSD_monitor(
    nx = 200,
    ny = 200,
    filename="pre_det",
    xwidth = 0.01,
    yheight = 0.01,
    restore_neutron = 1
) AT (0, 0, -thick) RELATIVE pre_samp_arm


COMPONENT arm_sample = Arm(
    
) AT (0, 0, 0.025) RELATIVE pre_samp_arm
ROTATED (0, 0, 0) RELATIVE pre_samp_arm


COMPONENT box_sample = Incoherent(
    xwidth=0.03,
    yheight=0.03,
    focus_xw=focus_w,
    focus_yh=focus_h,
    target_index=3,
    zdepth=thick,
    order=order
) WHEN (sample_case)
AT (0, 0, thick/2) RELATIVE arm_sample
ROTATED (0, 0, 0) RELATIVE arm_sample

COMPONENT total_scat = PSD_monitor_4PI(
    nx = 360,
    ny = 360,
    radius=0.1,
    restore_neutron=1,
    filename="total_scat"
) WHEN (total_scattering)
AT (0, 0, 0) RELATIVE arm_sample


COMPONENT arm_det = Arm(
) AT (0, 0, 0) RELATIVE arm_sample
ROTATED (0, det_rot, 0) RELATIVE arm_sample


COMPONENT psd_det = PSD_monitor(
    nx = 200,
    ny = 200,
    filename="post_det",
    xwidth = 0.031622776602,
    yheight = 0.031622776602,
    restore_neutron = 1
) AT (0, 0, 1) RELATIVE arm_det


/* This section is executed when the simulation ends (C code). Other    */
/* optional sections are : SAVE                                         */
FINALLY
%{
%}
/* The END token marks the instrument definition end */
END