/*******************************************************************************
*         McXtrace instrument definition URL=http://www.mcxtrace.org
*
* Instrument: NuSTAR_1shell
*
* %Identification
* Written by: Erik B Knudsen <erkn@fysik.dtu.dk> & Desiree D. M. Ferreira <desiree@space.dtu.dk> (email)
* Date: 12/12/2016
* Origin: DTU Physics/DTU Space
* Release: McXtrace 1.2
* Version: 1.0
* %INSTRUMENT_SITE: AstroX_NASA
*
* Single shell model of the NuSTAR-optic in use as telescope.
*
* %Description
* A model of the NuSTAR-telescope using just a single shell as optical element.
* That means to make use of this instrument
* it is necessary to run a series of simulation while varying the input parameter shellnumber.
*
* The model needs as input two geometry files (parabolic_datafile and hyperbolic_datafile), which contains (in ascii) tabled details about the geometry of the shells.
* The data in the geomfile is assumed to be in the format:
* #row  L/m  rad_h/m  rad_m/m  rad_p/m
* ...
* row: running index for the rows/rings
* L: The length of the plates for this ring
* rad_h: The radius at the "hyperbolic" end of the optic. At the detector end.
* rad_m: The radius at the midpoint of the optic. This is the reference.
* rad_p: The radius at the "parabolic" and of the optic. At the source end.
*
* Example: NuSTAR_1shell.instr shellnumber=1
*
* %Parameters
* FL: [m]            The focal length of the optical system
* FLd: [m]           The distance from the optics centre to the detector. If 0 the detector is placed at the focal plane.
* optics_dist: [m]   The distance between souce and optic. In space this would be quite large :-).
* offaxis_angle: [arcmin] Angle of collimated light from source
* reflectivity: [ ]  Data file containing reflectivities (such as from IMD)
* E0: [keV]          Central energy of X-rays
* dE: [keV]          Half spread of energy spectrum to be emitted from source
* shellnumber: [ ]   The row number for the miror module. This defines the shell.
* parabolic_datafile: [ ]  File which contains the geometry of the parabolic (primary) mirrors. (i.e. radii,lengths)
* hyperbolic_datafile: [ ] File which contains the geometry of the hyperbolic (secondary) mirrors. (i.e. radii,lengths)
* drx: [arcsec] Rotation of shell along X
* dry: [arcsec] Rotation of shell along Y
* drz: [arcsec] Rotation of shell along Z
*
* %Link
* <a href="http://www.cosmos.esa.int/web/athena">The ATHENA web pages @ ESA</a>
*
* %End
*******************************************************************************/

/* Change name of instrument and input parameters with default values */
DEFINE INSTRUMENT NuSTAR_1shell(FL=10.12,FLd=0, optics_dist=10,
        offaxis_angle=0, drx=0,dry=0,drz=0 ,
        string reflectivity="mirror_coating_unity.txt", E0=5, dE=0.001, int shellnumber=0,
        string parabolic_datafile="om_con_1a_110901_t1.txt", string hyperbolic_datafile="om_con_3a_110901_t1.txt")

/* The DECLARE section allows us to declare variables or  small      */
/* functions in C syntax. These may be used in the whole instrument. */
DECLARE
%{
    double RP,RMP,PL,RH,RMH,HL,PHP,PHH;
    int Pcoat,Hcoat;
    #pragma acc declare create(alphax)
    double alphax,alphay;

    double pore_width=0.83e-3;

    char *coatings[]={
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt",
        "mirror_coating_unity.txt"
    };

%}

USERVARS %{
  int parascatter;
  int hyperscatter;
  double pararef;
  double hyperref;
%}

/* The INITIALIZE section is executed when the simulation starts     */
/* (C code). You may use them as component parameter values.        */
INITIALIZE
%{
    if (offaxis_angle){
        alphax=offaxis_angle * MIN2RAD;
    }

    t_Table tpa,tph;
    int status=0;    
    if( (status=Table_Read(&tpa,parabolic_datafile,0))<1){
        fprintf(stderr,"Error: Could not open/parse file %s\n",parabolic_datafile);
        exit(-1);
    } 
    if( (status=Table_Read(&tph,hyperbolic_datafile,0))<1){
        fprintf(stderr,"Error: Could not open/parse file %s\n",hyperbolic_datafile);
        exit(-1);
    } 

    /*extract relevant numbers from the tables*/
    RP=Table_Index(tpa,shellnumber,0)*1e-3;
    RMP=Table_Index(tpa,shellnumber,1)*1e-3;
    PL=Table_Index(tpa,shellnumber,5)*1e-3;
    Pcoat=Table_Index(tpa,shellnumber,7);
    RH=Table_Index(tph,shellnumber,1)*1e-3;
    RMH=Table_Index(tph,shellnumber,0)*1e-3;
    HL=Table_Index(tph,shellnumber,5)*1e-3;
    Hcoat=Table_Index(tph,shellnumber,7);
    PHH=(Table_Index(tph,shellnumber,9)-Table_Index(tph,shellnumber,8))*1e-3;
    PHP=(Table_Index(tpa,shellnumber,9)-Table_Index(tpa,shellnumber,8))*1e-3;

    if(!FLd){
        FLd=FL;
    }
    #pragma acc update device(alphax)
%}
/* instrument is defined as a sequence of components.  */
TRACE

/* The Arm() class component defines reference points and orientations  */
/* in 3D space. Every component instance must have a unique name. Here, */
/* Origin is used. This Arm() component is set to define the origin of  */
/* our global coordinate system (AT (0,0,0) ABSOLUTE). It may be used   */
/* for further RELATIVE reference, Other useful keywords are : ROTATED  */
/* EXTEND GROUP PREVIOUS. Also think about adding an xray source !    */
/* Progress_bar is an Arm displaying simulation progress.               */
COMPONENT Origin = Progress_bar()
AT (0,0,0) ABSOLUTE
EXTEND
%{
    parascatter=0;
    hyperscatter=0;
%}

COMPONENT src = Source_div(
        xwidth=0,yheight=2.0*PHP,focus_aw=0,focus_ah=0,E0=E0,dE=dE)
AT(0,RP-PHP/2.0,0) RELATIVE Origin

COMPONENT srcradsym = Arm()
AT(0,0,0) RELATIVE Origin
EXTEND
%{
    do {
        x=0;
        double eta=2*M_PI*rand01();
        rotate(x,y,z, x,y,z, eta, 0,0,1);
    } while(0);
%}


COMPONENT srcoffaxis= Arm()
WHEN(alphax!=0) AT(0,0,0) RELATIVE Origin
ROTATED (0,0,0) RELATIVE Origin
EXTEND
%{
    do {    
        rotate(kx,ky,kz, kx,ky,kz, alphax, 0,1,0);
        x-=INSTRUMENT_GETPAR(optics_dist)*sin(alphax);
        SCATTER;
    }while(0);
%}

COMPONENT detector_pre_optics = PSD_monitor(restore_xray=1, xwidth=3, yheight=1.5, nx=101, ny=51, filename="det_preo.dat")
AT(0,0,optics_dist) RELATIVE Origin

COMPONENT optics_centre = Arm()
AT(0,0,optics_dist) RELATIVE Origin

COMPONENT misalign = Arm()
AT(0,0,0) RELATIVE optics_centre
ROTATED (drx/3600.0, dry/3600.0, drz/3600.0) RELATIVE optics_centre

COMPONENT Shell_p_1 = Shell_p(
        radius_p=RP, radius_m=RMP, zdepth=PL, Z0=FL, yheight=PHP, mirror_reflec=coatings[Pcoat], R_d=0)
AT(0,0,-1e-4) RELATIVE misalign
EXTEND
%{
    if (SCATTERED){
        parascatter=SCATTERED;
    }
%}

COMPONENT midopdet = PSD_monitor(
        restore_xray=1,xwidth=1,yheight=1,nx=201,ny=101, filename="midop.dat")
AT(0,0,0) RELATIVE misalign

COMPONENT Shell_h_1 = Shell_h(
        radius_m=RMH, radius_h=RH, zdepth=HL, Z0=FL, yheight=PHH, mirror_reflec=coatings[Hcoat], R_d=0)
AT(0,0,0) RELATIVE misalign
GROUP hyperoptics
EXTEND
%{
    if (SCATTERED){
        hyperscatter=SCATTERED;
    }
%}

COMPONENT big_detector = PSD_monitor(restore_xray=1, xwidth=0.2, yheight=0.2, nx=201, ny=201, filename="bigdet.dat")
AT(0,0,FLd) RELATIVE optics_centre

COMPONENT big_detector_parab = copy(big_detector)(filename="bigdet_parab.dat")
WHEN(parascatter>1 && hyperscatter<=1) AT(0,0,FLd) RELATIVE optics_centre

COMPONENT big_detector_hyper = COPY(big_detector)(filename="bigdet_hyperb.dat")
WHEN(hyperscatter>1 && parascatter<=1) AT(0,0,FLd) RELATIVE optics_centre

COMPONENT focal_detector = PSD_monitor(restore_xray=1,xwidth=1e-2, yheight=1e-2, nx=201, ny=201, filename="focal_det.dat")
AT(0,0,FLd) RELATIVE optics_centre
COMPONENT superfocal_detector = PSD_monitor(restore_xray=1,xwidth=1e-6, yheight=1e-6, nx=201, ny=201, filename="superfocal_det.dat")
AT(0,0,FLd) RELATIVE optics_centre
COMPONENT ultrafocal_detector = PSD_monitor(restore_xray=1,xwidth=1e-12, yheight=1e-12, nx=201, ny=201, filename="ultrafocal_det.dat")
AT(0,0,FLd) RELATIVE optics_centre

/* 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