/*******************************************************************************
*         McXtrace instrument definition URL=http://www.mcxtrace.org
*
* Instrument: ATHENA_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_ESA
*
* Single shell model of the ATHENA SPO-optic in use as telescope.
*
* %Description
* A model of the ATHENA-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 a file geomfile, 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  width/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.
* width: pore width ?
*
* Example: ATHENA_1shell.instr shellnumber=1
*
* %Parameters
* FL: [m] The focal length of the optical system
* optics_dist: [m] The distance between souce and optic. In space this would be quite large :-).
* SRC_POS_X: [m] Displacement of source along X
* SRC_POS_Y: [m] Displacement of source along Y
* 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.
* geomfile: [ ] File which contains the geometry of the pores (i.e. radii,lengths)
*
* %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 ATHENA_1shell(FL=12, optics_dist=10,
        SRC_POS_X=0, SRC_POS_Y=0, offaxis_angle=0, drx=0,dry=0,drz=0 ,
        string reflectivity="mirror_coating_unity.txt", E0=5, dE=0.001, int shellnumber=1,
	string geomfile="ATHENA_rings_1_20.dat")
DEPENDENCY "-DDATPATH=GETPATH(data/)"

/* The DECLARE section allows us to declare variables or  small      */
/* functions in C syntax. These may be used in the whole instrument. */
DECLARE
%{
    double PX[1024],PY[1024], PR[1024], PA[1024];
    int IDX[2014],ROW[2014], RIDX[1024];
    double PRP[1024],PRM[1024],PRH[1024],PRL[1024],PRW[1024];

    double src_pos_x;
    double src_pos_y;
    #pragma acc declare create(alphax)
    double alphax,alphay;

    double pore_width=0.83e-3;
    double pore_height=0.605e-3;

    #define QUOTE(name) #name
    #define STR(macro) QUOTE(macro)

    #ifndef DATPATH
    #define DATPATH=/usr/share/mcxtrace/3.1/data/
    #endif

    #define DATAPATH STR(DATPATH)

    char fullpath_geomfile[128];
%}

USERVARS %{
  double hyperref;
  double pararef;
  int parascatter;
  double hyperscatter;
  double pstore;
  long long nid;
%}

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

    sprintf(fullpath_geomfile,"%s/%s",DATAPATH,geomfile);

    FILE *fp=fopen(fullpath_geomfile,"rb");
    if (fp==NULL){
        fprintf(stderr,"Error: Cannot open file \'%s\'. Aborting.\n",fullpath_geomfile);
        exit(-1);
    }
    int ii=0;
    char line[1024];
    fgets(line,1023,fp);
    ii=0;
    while (!feof(fp)){
        fscanf(fp,"%d %lf %lf %lf %lf %lf\n",RIDX+ii,PRL+ii, PRH+ii, PRM+ii, PRP+ii, PRW+ii);
        ii++;
    }
    fclose(fp);
    #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;
    nid++;
    hyperref=1;
    pararef=1;
%}

COMPONENT src = Source_div(
        xwidth=0,yheight=2.0*pore_height,focus_aw=0,focus_ah=0,E0=E0,dE=dE)
AT(0,PRP[shellnumber-1]-pore_height/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
EXTEND
%{
    pstore=p;
%}

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=PRP[shellnumber-1], radius_m=PRM[shellnumber-1], zdepth=PRL[shellnumber-1], Z0=FL, yheight=0.605e-3, mirror_reflec=reflectivity, R_d=0)
AT(0,0,0) RELATIVE misalign
EXTEND
%{
    if (SCATTERED){
        parascatter=SCATTERED;
        pararef=p/pstore;
    }
%}

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

COMPONENT Shell_h_1 = Shell_h(
    radius_m=PRM[shellnumber-1], radius_h=PRH[shellnumber-1], zdepth=PRL[shellnumber-1], Z0=FL, yheight=0.605e-3, mirror_reflec=reflectivity, R_d=0)
AT(0,0,0) RELATIVE misalign
GROUP hyperoptics
EXTEND
%{
    if (SCATTERED){
        hyperscatter=SCATTERED;
        hyperref=p/(pstore*pararef);
    }
%}

COMPONENT fourpi = PSD_monitor_4PI(restore_xray=1,filename="sphere.dat", radius=2, nx=51,ny=51)
AT(0,0,0) RELATIVE optics_centre 

COMPONENT detector_post_optics = PSD_monitor(restore_xray=1,xwidth=1.5, yheight=1, nx=201, ny=101, filename="det_posto.dat")
AT(0,0,optics_dist+0.5) RELATIVE Origin

COMPONENT paramon = Monitor_nD(
        restore_xray=1,filename="paramond",xwidth=3, yheight=1.5, options="list=2000 user1,user2,user3", user1="nid",user2="parascatter",user3="pararef",
        username1="ray_id",username2="parabolic_shell_id",username3="parabolic_reflectivity")
AT(0,0,optics_dist+0.5) RELATIVE Origin

COMPONENT hypermon = Monitor_nD(
        restore_xray=1,filename="hypermond",xwidth=3, yheight=1.5, options="list=2000 user1,user2,user3", user1="nid",user2="hyperscatter",user3="hyperref",
        username1="ray_id",username2="hyperbolic_shell_id",username3="hyperbolic_reflectivity")
AT(0,0,optics_dist+0.5) RELATIVE Origin

COMPONENT big_detector = PSD_monitor(restore_xray=1, xwidth=0.16, yheight=1, nx=201, ny=201, filename="bigdet.dat")
AT(0,0,FL) 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,FL) 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,FL) 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,FL) RELATIVE optics_centre

COMPONENT FLmond= Monitor_nD(
        restore_xray=1,filename="FLmond",xwidth=0.1, yheight=.1, options="list=all user1 x y k E", user1="nid",
        username1="ray_id")
AT(0,0,FL) 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