/*******************************************************************************
 * 	McStas simulation of the Monochromator Chamber of Necsa's Neutron Diffraction Facility
 *
 * Instrument: SAFARI-1 Beam Ports
 *
 * %Identification
 * Written by: Deon Marais (deon.marais@necsa.co.za)
 * Adapted by: Tjatji Tjebane (tjatji.tjebane@necsa.co.za)
 * Date: June 2021
 * Origin: Necsa
 * Release: McStas 2.7
 * Version: $Revision 0.1 $
 * %INSTRUMENT_SITE: Necsa
 *
 * SAFARI-1 Beam Ports
 *
 * %Description
 * The Necsa SAFARI-1 research reactor, South Africa. Beam ports 1(SANS) 2(NRAD) 5(NDIFF)
 * 
 * %Parameters
 * Input parameters		Unit		Description
 * ----------------		----		-----------
 * source_lam_min:		[Angs] 		Minimum wavelenth of source
 * source_lam_max:		[Angs]		Maximum wavelenth of source
 * fluxFile:			[string]	Flux file location
 *
 * %Link
 * The South African Nuclear Energy Corporation <a href="http://www.necsa.co.za"> website</a>	
 *
 * %End
 *******************************************************************************/

DEFINE INSTRUMENT BEAM_PORTS(source_lam_min=0.1, source_lam_max=30.0, NVS_torsion_angle =12.03211, NVS_nu=477.4648293,
				char *fluxFile="SAFARI-BP1.mcpl", MCNP_Normalisation = 1.4902506E+18, MCNP_Histories = 57600759, int repetition = 1
			) 

DECLARE
%{
char *fluxFile;
double MCNP_norm;
double MCNP_hist;
double source_e_min;
double source_e_max;

char op_PSDcm2Sq[256];
char op_PSDcm2Cr[256];
char op_PSD[256];
char op_FLUX[256];
char op_LAM[256];
char op_E[256];


double BP1Diam;

double Centre_to_BShield;
double B1Vessel_to_BShield;
double Water_width;
double PShutterCavityWidth;
double PShutterCavityLength;
double PShutterCavityHeight;

double Be_filter_length;
double Filter_cavity_length;
double NVS_length;
//double NVS_nu;
double NVS_torsion_angle;
double DetRadius;
double Collimator_length;
double Collimator_diameter;

%}

INITIALIZE
%{
MCNP_norm = MCNP_Normalisation;
MCNP_hist = MCNP_Histories;

BP1Diam = 17.46/100.0;	//cm to m

Centre_to_BShield = (4430.0+(828.68-415.05))/1000.0;  //mm to m
B1Vessel_to_BShield = 4.0431; //4,8768 - 0.8337
Water_width = 1.180;
PShutterCavityWidth = 78.4/100.0;
PShutterCavityLength = 77.0/100.0;
PShutterCavityHeight = 76.5/100.0;

Be_filter_length = 20.0/100.0;
Filter_cavity_length = 76.0/100.0;
NVS_length = 64.8/100.0;
//NVS_nu = 477.4648293;
//NVS_torsion_angle = 12.03211;
DetRadius = 0.5;
Collimator_length = 5.0;
Collimator_diameter = 0.05;

source_e_min = (9.044888059/source_lam_max)*(9.044888059/source_lam_max);	
source_e_max = (9.044888059/source_lam_min)*(9.044888059/source_lam_min);	

//source_e_min = (0.286/source_lam_max)*(0.286/source_lam_max)/1e4;		
//source_e_max = (0.286/source_lam_min)*(0.286/source_lam_min)/1e4;

printf ("\n------------------\n");
printf ("\n------------------\nENERGY RANGES: %g to %g\n",source_e_min,source_e_max);

sprintf(op_PSDcm2Sq, "parallel, verbose, square, x limits=[-0.1 0.1] bins=100, y limits=[-0.1 0.1] bins = 100, cm2");
sprintf(op_PSDcm2Cr, "parallel, verbose, disk, x limits=[-0.1 0.1] bins=100, y limits=[-0.1 0.1] bins = 100, cm2");
sprintf(op_PSD, "parallel, verbose, square, x limits=[-0.1 0.1] bins=100, y limits=[-0.1 0.1] bins = 100");
sprintf(op_E, "parallel, verbose, square, energy limits==[%g %g] bins=100, cm2", 0.0, 300.0);
sprintf(op_LAM, "parallel, verbose, square, lambda limits=[%g %g] bins=100, cm2", source_lam_min, source_lam_max);



%}

TRACE
//**************************************************************************************************
//Source
COMPONENT Progress = Progress_bar(percent=1,flag_save=0) 
	AT (0, 0, 0) ABSOLUTE
	
COMPONENT CoreCenter = Arm()
 AT (0, 0, 0) ABSOLUTE

COMPONENT MCNP_Axis = Arm()
 AT  (0, 0, 0) RELATIVE CoreCenter
 ROTATED (-90, 0, 180) ABSOLUTE

COMPONENT B1Angle = Arm()
  AT (0, 0.1143, 0) ABSOLUTE
  ROTATED (0, 74.34, 0) ABSOLUTE	//From CAD model
COMPONENT B1Vessel = Arm()
  AT (0, 0, 0.8337) Relative B1Angle

/*COMPONENT B2Angle = Arm()
  AT (0, -0.14922, -0.15243) ABSOLUTE
  ROTATED (0, 51.71, 0) ABSOLUTE	//From CAD model
COMPONENT B2Vessel = Arm()
  AT (0, 0, 0.9233) Relative B2Angle

COMPONENT B5Angle = Arm()
  AT (0, 0.1143, 0) ABSOLUTE
  ROTATED (0, -10.62, 0) ABSOLUTE	//From CAD model
COMPONENT B5Vessel = Arm()
  AT (0.0, 0.0, 0.8337) Relative B5Angle
 */
//MCPL input
COMPONENT MCNPSource = MCPL_input(filename = fluxFile,E_smear=0.00000, pos_smear=0.0000, dir_smear=0.000, verbose=1, repeat_count=repetition)
AT (0,0,0) RELATIVE MCNP_Axis
EXTEND %{
  p*=MCNP_norm;
  p/=MCNP_hist;
%}



//----------------------------------------------------------------------------
// Beam 1

//COMPONENT B2Stop = Beamstop (radius = BP2Diam/2.0)
//  AT (0, 0, 0.02) RELATIVE B2Vessel

//COMPONENT B5Stop = Beamstop (radius = BP5Diam/2.0)
//  AT (0, 0, 0.02) RELATIVE B5Vessel

/*COMPONENT PSDcm2Cr_Vessel= Monitor_nD(
    filename = "PSDcm2Cr_Vessel.sim", restore_neutron = 1,
      xwidth = BP1Diam,
      options = op_PSDcm2Cr)
  AT (0, 0, 0.02) RELATIVE B1Vessel 
*/  
  COMPONENT PSD_Vessel= Monitor_nD(
      filename = "PSD_Vessel.sim", restore_neutron = 1,
        xwidth = 0.1, yheight = 0.1,
        options = op_PSD)
  AT (0, 0, 0.02) RELATIVE B1Vessel 

COMPONENT LAM_Vessel = Monitor_nD(
      filename = "LAM_Vessel.sim", restore_neutron = 1, options = op_LAM,
      xwidth = 0.1, yheight = 0.1)
  AT (0, 0, 0.02) RELATIVE B1Vessel

 
//Quartz Filter
COMPONENT SiO2_filter = PowderN(reflections="SiO2_quartza.laz", radius = 0.01, thickness = 0.001,
concentric = 1, p_interact=0.1)
AT (0,0,Water_width +0.1) RELATIVE B1Vessel


COMPONENT SiO2_PSDmonitor=PSD_monitor( xwidth = BP1Diam/2.0, yheight = BP1Diam/2.0, nx=100, ny=100, filename="PSD_SiO2_filter.psd")
  AT (0, 0, 0.01) RELATIVE SiO2_filter 

COMPONENT LAM_SiO2_filter = Monitor_nD(
      filename = "LAM_SiO2_filter.sim", restore_neutron = 1, options = op_LAM,
      xwidth = 0.1, yheight = 0.1)
  AT (0, 0, 0.02) RELATIVE SiO2_filter 


/*
//Bismuth Filter
COMPONENT Bi_filter = PowderN(reflections="Bi.laz", radius = 0.01, thickness = 0.001,
concentric = 1, p_interact=0.1)
AT (0,0,B1Vessel_to_BShield+0.1) RELATIVE B1Vessel

//COMPONENT Bi_PSDmonitor=PSD_monitor( xwidth = BP1Diam, yheight = BP1Diam, nx=100, ny=100, filename="PSD_Bi_filter.psd")
//  AT (0, 0, 0.01) RELATIVE Bi_filter


COMPONENT PSDcm2Cr_Bi_filter= Monitor_nD(
    filename = "PSDcm2Cr_Bismuth_filter.sim", restore_neutron = 1,
      xwidth = 0.06,
      options = op_PSDcm2Cr)
  AT (0, 0, Filter_cavity_length+0.01) RELATIVE Bi_filter 


COMPONENT LAM_Bi_filter = Monitor_nD(
      filename = "LAM_Bismuth_filter.sim", restore_neutron = 1, options = op_LAM,
      xwidth = 0.1, yheight = 0.02)
  AT (0, 0, 0.02) RELATIVE Bi_filter
*/

//Beryllium Filter

COMPONENT Berylium_filter = Filter_gen(xwidth=0.06, yheight=0.06, thickness=Be_filter_length/0.05, filename="Be.trm")
	AT (0, 0, B1Vessel_to_BShield+0.4) RELATIVE  B1Vessel 

/*COMPONENT Be_filter = Filter_gen(
  filename = "Be.trm", options = "multiply",
  xwidth=0.06, yheight=0.06, thickness=4) 
  AT (0, 0, B1Vessel_to_BShield+0.4) RELATIVE B1Vessel
*/


//COMPONENT Be_PSDmonitor=PSD_monitor( xwidth = BP1Diam, yheight = BP1Diam, nx=100, ny=100, filename="PSD_Be_filter.psd")
//  AT (0, 0, 0.01) RELATIVE Be_filter

  
COMPONENT PSDcm2Cr_Be_filter= Monitor_nD(
    filename = "PSDcm2Cr_Berylium_filter.sim", restore_neutron = 1,
      xwidth = 0.06,
      options = op_PSDcm2Cr)
  AT (0, 0, Filter_cavity_length+0.01) RELATIVE Berylium_filter 


//COMPONENT LAM_Berylium_filter= L_monitor(xmin=-0.1, xmax=0.1, ymin=-0.1, ymax=0.1, nL=20, filename="Be_Output.L", Lmin=2, Lmax=10)
//AT (0, 0, Filter_cavity_length+0.01) RELATIVE Berylium_filter

COMPONENT LAM_Berylium_filter = Monitor_nD(
      filename = "LAM_Be_filter.sim", restore_neutron = 1, options = op_LAM,
      xwidth = 0.1, yheight = 0.1)
  AT (0, 0, 0.02) RELATIVE Berylium_filter




//Velocity Selector
COMPONENT NVS = Selector( xwidth = 0.08, yheight = 0.06,
	length=NVS_length, nslit=70, d=0.0005, radius=0.21, alpha=NVS_torsion_angle, 
	nu=NVS_nu)
  AT (0, 0, B1Vessel_to_BShield+Filter_cavity_length+0.1) RELATIVE B1Vessel


/*
// ILL NVS

COMPONENT NVS = Selector(xmin=-0.015,  xmax=0.015,  ymin=-0.025, ymax=0.025, length=0.25,
	nslit=72,d=0.0004, radius=0.12, alpha=48.298, nu=500)
  AT (0, 0, B1Vessel_to_BShield+0.1) RELATIVE B1Vessel
*/

/*COMPONENT NVS_PSDmonitor=PSD_monitor( xwidth = BP1Diam/2.0, yheight = BP1Diam/2.0, nx=100, ny=100, filename="PSD_NVS.psd")
  AT (0, 0, 0.01) RELATIVE NVS
 */

COMPONENT PSDcm2Cr_NVS= Monitor_nD(
    filename = "PSDcm2Cr_NVS.sim", restore_neutron = 1,
      xwidth = 0.08,
      options = op_PSDcm2Cr)
  AT (0, 0, 0.01) RELATIVE NVS 


COMPONENT LAM_NVS = Monitor_nD(
      filename = "LAM_NVS.sim", restore_neutron = 1, options = op_LAM,
      xwidth = 0.1, yheight = 0.1)
  AT (0, 0, 0.02) RELATIVE NVS 
  
  
// Collimator  
COMPONENT SampleCollimator = Collimator_linear(xmin=-Collimator_diameter/2.0, xmax=Collimator_diameter/2.0, ymin=-Collimator_diameter/2.0, 
	ymax=Collimator_diameter/2.0, length=Collimator_length, divergence=0)
 AT (0, 0, B1Vessel_to_BShield+Filter_cavity_length+NVS_length+0.1) RELATIVE B1Vessel

//COMPONENT Collimator_PSDmonitor=PSD_monitor( xwidth = BP1Diam/2.0, yheight = BP1Diam/2.0, nx=100, ny=100, filename="PSD_Collimator.psd")
//  AT (0, 0, 0.01) RELATIVE  SampleCollimator


/*COMPONENT PSD_SampleCollimator= Monitor_nD(
      filename = "PSD_SampleCollimator.sim", restore_neutron = 1, options = op_PSD,
        xwidth = 0.2, yheight = 0.2)
  AT (0, 0, 0.02) RELATIVE SampleCollimator 
*/

COMPONENT PSDcm2Cr_SampleCollimator= Monitor_nD(
    filename = "PSDcm2Cr_SampleCollimator.sim", restore_neutron = 1,
      xwidth = 0.08,
      options = op_PSDcm2Cr)
  AT (0, 0, 0.01) RELATIVE SampleCollimator 

COMPONENT LAM_SampleCollimator = Monitor_nD(
      filename = "LAM_SampleCollimator.sim", restore_neutron = 1, options = op_LAM,
      xwidth = 0.1, yheight = 0.1)
  AT (0, 0, 0.02) RELATIVE SampleCollimator 

 
  


FINALLY
%{
%}
END