/*******************************************************************************
*         McStas instrument definition URL=http://www.mcstas.org
*
* Instrument: BTsimple.instr
*
* %I
* Written by: Peter Willendrup <pkwi@fysik.dtu.dk>
* Date: 2018/03/21
* Origin: DTU
* %INSTRUMENT_SITE: Templates
*
* Example instrument showing how to calculate brilliance transfer using L_monitor's and WHEN statements
*
* %D
* Example instrument showing how to calculate brilliance transfer using L_monitor's and WHEN statements
*
* %Example: -n1e7 BTsimple.instr lambda=10 Detector: BTransfer_I=75.4758
*
* %P
* lambda:  [AA]  Central wavelength produced at the source
* dlambda: [AA]  Halfwidth of wavelength band produced at the source
* maxhd:   [deg] Maxiamal horzizontal divergence accepted
* maxvd:   [deg] Maxiamal vertical divergence accepted
* gw:      [m]   Opening-width of ellipitic guide
* gh:      [m]   Opening-height of ellipitic guide
* gL:      [m]   Length of ellipitic guide
* gm:      [1]   M-value of the guide coating
* delta1:  [m]   Optional focal displacement of guide wrt source
* delta2:  [m]   Optional focal displacement of guide wrt sample position
* Nbins:   [1]   Number of bins on Wavelength / Brilliance transfer monitors
*
* %L
*
* %E
*******************************************************************************/
DEFINE INSTRUMENT BTsimple(lambda=10, dlambda=9.9, maxhd=3, maxvd=3, gw=0.1, gh=0.1, gL=50, gm=6, delta1=0, delta2=0, int Nbins=101)

DECLARE
%{	
  /* Arrays for storing the calculated brilliance transfer */
  DArray1d BT_N;
  DArray1d BT_p;
  DArray1d BT_p2;
%}

USERVARS
%{
  /* Variables for calculation of neutron divergence  */
  double VertDiv;
  double HorDiv;
%}

INITIALIZE
%{
  BT_N = create_darr1d(Nbins);
  BT_p = create_darr1d(Nbins);
  BT_p2 = create_darr1d(Nbins);
%}

TRACE

COMPONENT origin = Progress_bar()
AT (0, 0, 0) RELATIVE ABSOLUTE
EXTEND %{
  /* Initialize divergencies to 0 */
  VertDiv = 0;
  HorDiv  = 0;
%}

COMPONENT source = Source_gen(
    /* ~ ILL cold source */
    T1=413.5,I1=10.22e12,T2=145.8,I2=3.44e13,
    T3=40.1 ,I3=2.78e13,
    dist=2,
    radius=0.06,
    focus_xw=gw, 
    focus_yh=gh, 
    lambda0=lambda, 
    dlambda=dlambda)
AT (0, 0, 0) RELATIVE PREVIOUS

COMPONENT CalcDiv1 = Arm()
AT (0,0,0) RELATIVE PREVIOUS
EXTEND %{
  /* Calculate divergence of ray as emitted by the source */
  VertDiv = RAD2DEG*atan2(vy,vz);
  HorDiv  = RAD2DEG*atan2(vx,vz);
%}

/* Measure incoming phase-space */
COMPONENT BT_in = L_monitor(xwidth=0.02, yheight=0.02, filename="BT_in.dat",Lmin=lambda-dlambda, Lmax=lambda+dlambda, nL=Nbins, restore_neutron=1)
WHEN ((VertDiv <= maxvd) && (HorDiv <=  maxhd)) AT (0,0,2.0) RELATIVE  origin


/* Visualise incoming divergence */
COMPONENT Div_in = Divergence_monitor(nh=20, nv=20, filename="Div_in",
				      xwidth=gw,yheight=gh, maxdiv_h=5, maxdiv_v=5, restore_neutron=1)
AT (0,0,2.0) RELATIVE  origin

/* Visualise incoming beam footprint */
COMPONENT PSD_in = PSD_monitor(ny=20, nx=20, filename="PSD_in",
				      xwidth=gw,yheight=gh,restore_neutron=1)
AT (0,0,2.0) RELATIVE  origin

/* Our elliptical guide */
  COMPONENT guide = Elliptic_guide_gravity(
    l=gL,
    linxw = 2+delta1,
    loutxw= 2+delta2,
    linyh = 2+delta1,
    loutyh= 2+delta2,
    xwidth=gw,
    yheight=gh,
    alpha=0,
    m=gm,
    W=0)
AT (0, 0, 2) RELATIVE origin

COMPONENT CalcDiv2 = Arm()
AT (0,0,gL) RELATIVE PREVIOUS
EXTEND %{
    /* Calculate divergence of ray after the guide */
    VertDiv = RAD2DEG*atan2(vy,vz);
    HorDiv = RAD2DEG*atan2(vx,vz);
%}

/* Measure outgoing phase-space @ sample position */
COMPONENT BT_out = L_monitor(xwidth=0.02, yheight=0.02, filename="BT_out.dat",Lmin=lambda-dlambda, Lmax=lambda+dlambda, nL=Nbins, restore_neutron=1)
WHEN ((VertDiv <= maxvd) && (HorDiv <=  maxhd)) AT (0,0,2.0+gL) RELATIVE  guide


/* Visualise outgoing divergence */
COMPONENT Div_out = Divergence_monitor(nh=20, nv=20, filename="Div_out",
				      xwidth=gw,yheight=gh, maxdiv_h=5, maxdiv_v=5, restore_neutron=1)
AT (0,0,2.0+gL) RELATIVE  guide

/* Visualise outgoing beam footprint */
COMPONENT PSD_out = PSD_monitor(ny=20, nx=20, filename="PSD_out",
				      xwidth=gw,yheight=gh,restore_neutron=1)
AT (0,0,0) RELATIVE  PREVIOUS

FINALLY
%{
/* This adds another "monitor" that measures BT_out / BT_in */
/* In MPI-mode do the calculation on master only.           */
#ifdef USE_MPI
  if(mpi_node_rank == mpi_node_root) {
#endif
    int j;
    double* tmpN;
    double* tmpp1;
    double* tmpp2;
    double* tmpd1;
    double* tmpd2;
    tmpN =COMP_GETPAR(BT_out,L_N);
    tmpp1=COMP_GETPAR(BT_in, L_p);
    tmpp2=COMP_GETPAR(BT_out,L_p);
    tmpd1=COMP_GETPAR(BT_in, L_p2);
    tmpd2=COMP_GETPAR(BT_out,L_p2);
    
    for (j=0;j<Nbins;j++) {
      BT_N[j]=tmpN[j];
      if (tmpp1[j] != 0) {
	BT_p[j]=tmpp2[j]/tmpp1[j];
      } else {
	BT_p[j]=0;
      }
      if ((tmpp1[j] != 0) && (tmpp2[j] != 0)) {
	BT_p2[j]=sqrt((tmpd1[j]/tmpp1[j])*(tmpd1[j]/tmpp1[j])  + (tmpd2[j]/tmpp2[j])*(tmpd2[j]/tmpp2[j]));
      } else
	BT_p2[j]=0;
    }
#ifdef USE_MPI
  }
  // Non-master MPI nodes wait here before we write the file to disk
  MPI_Barrier(MPI_COMM_WORLD);
#endif
  // we can not call the DETECTOR_OUT_1D as it transparently makes use of
  // NAME_CURRENT_COMP and POS_A_CURRENT_COMP, which are set to match '_comp'
  Rotation Rot;
  rot_set_rotation(Rot,0,0,0);
  mcdetector_out_1D(
		    "Brilliance transfer",
		    "Wavelength [AA]",
		    "BT",
		    "L", lambda-dlambda, lambda+dlambda, Nbins,
		    &BT_N[0],&BT_p[0],&BT_p2[0],
		    "Brilliance_transfer", "BTransfer", coords_set(0,0,0),Rot,9999);

%}

END