/*******************************************************************************
*         McXtrace instrument definition URL=http://www.mcxtrace.org
*
* Instrument: ATHENA_1mm
*
* %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 pore version of the ATHENA SPO-optic in use as telescope.
*
* %Description
* A model of the ATHENA-telescope using just a single mirror module 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 porenumber.
* At present, the porenumber really means mirror module number.
*
* The model needs as input two files positionfile and geomfile, which contain (in ascii) tabled details about the positions of mirror
* modules and their geometry.
* The data in the position file is assumed to be columnar with the following format:
* #index   row    X/m    Y/m    R/m    alpha/deg.
* ...
*
* index: a running number identifying the mirror module within a petal.
* row: the a running number for which row (ring) this mirror module belongs to (n.b. this is used as an index inot the geomfile).
* X: The x-coordinate of the mirror module (ignored because redundant with R,alpha)
* Y: The y-coordinate of the mirror module (ignored because redundant with R,alpha)
* R: The radial position of the mirror module
* alpha: The angular position of the mirror module
*
* Similarly the format of the geomfile is:
* #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_1mm.instr porenumber=3
*
* %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
* porenumber: [ ] Actually the mirror module number.
* positionfile: [ ] File containing the pore/mirror module positions.
* geomfile: [ ] File which contains the geometry of the pores (i.e. radii,lengths)
* dPx:      [m] Offset/displacement of parabolic pore along x from its theoretical position.
* dPy:      [m] Offset/displacement of parabolic pore along y from its theoretical position.
* dPz:      [m] Offset/displacement of parabolic pore along z from its theoretical position.
* dPrx:     [arcsec] Rotational misalignment of parabolic pore around x.
* dPry:     [arcsec] Rotational misalignment of parabolic pore around y.
* dPrz:     [arcsec] Rotational misalignment of parabolic pore around z.
* dHx:      [m] Offset/displacement of hyperbolic pore along x from its theoretical position.
* dHy:      [m] Offset/displacement of hyperbolic pore along y from its theoretical position.
* dHz:      [m] Offset/displacement of hyperbolic pore along z from its theoretical position.
* dHrx:     [arcsec] Rotational misalignment of hyperbolic pore around x.
* dHry:     [arcsec] Rotational misalignment of hyperbolic pore around y.
* dHrz:     [arcsec] Rotational misalignment of hyperbolic pore around z.
* Hyper:    [ ] If non-zero the secondary mirror (hyperbolic) is active. Useful for debugging.
* Para:     [ ] If non-zero the primary mirror is acive (parabolic) is active. Useful for debugging.
* ringfile: [ ] File which contains deatiled plate descriptions.
*
* %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_1mm(FL=12, optics_dist=10,
        SRC_POS_X=0, SRC_POS_Y=0, offaxis_angle=0, 
        dPx=0,dPy=0,dPz=0, dPrx=0, dPry=0, dPrz=0,
        dHx=0,dHy=0,dHz=0, dHrx=0, dHry=0, dHrz=0,
        string reflectivity="mirror_coating_unity.txt", 
        E0=5, dE=0.001, int porenumber=1,string ringfile="ref_design_breaks.txt", 
        string positionfile="ATHENA_mm_1petal.dat", string geomfile="ATHENA_rings_1_20.dat", int Hyper=1, int Para=1)
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_positionfile[128];
    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_positionfile,"%s/%s",DATAPATH,positionfile);

    FILE *fp=fopen(fullpath_positionfile,"rb");
    if (fp==NULL){
        fprintf(stderr,"Error: Cannot open file \'%s\'. Aborting.\n",fullpath_positionfile);
        exit(-1);
    }
    int ii=0;
    char line[1024];
    fgets(line,1023,fp);
    while (!feof(fp)){
        fscanf(fp,"%d %d %lf %lf %lf %lf\n",IDX+ii, ROW+ii, PX+ii, PY+ii, PR+ii, PA+ii);
        ii++;
    }
    fclose(fp);

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

    fp=fopen(fullpath_geomfile,"rb");
    if (fp==NULL){
        fprintf(stderr,"Error: Cannot open file \'%s\'. Aborting.\n",fullpath_geomfile);
        exit(-1);
    }
    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 srca = Arm()
AT(0,0,0) RELATIVE Origin
ROTATED(0,0,-90+PA[porenumber-1]) RELATIVE Origin

COMPONENT src = Source_div(
        xwidth=100.0*pore_width,yheight=100.0*pore_height,focus_aw=0,focus_ah=0,E0=E0,dE=dE)
AT(0,PR[porenumber-1],0) RELATIVE srca

COMPONENT srcoffaxis= Arm()
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=3, nx=51, 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 a_1 = Arm()
AT(0,0,0) RELATIVE optics_centre
ROTATED (0,0,-90+PA[porenumber-1]) RELATIVE optics_centre

COMPONENT misalign_rot_p = Arm()
AT(dPx,dPy+PRM[ROW[porenumber-1]-1],dPz) RELATIVE a_1
ROTATED (dPrx/3600.0,dPry/3600.0,dPrz/3600.0) RELATIVE a_1

COMPONENT ref_p = Arm()
AT( 0,-PRM[ROW[porenumber-1]-1],0) RELATIVE misalign_rot_p

/*GROUP paraoptics*/

COMPONENT mm_p_1 = MM_p(
        pore_th=0, ring_nr=ROW[porenumber-1], Z0=FL, pore_width=pore_width , pore_height=pore_height, mirror_reflec=reflectivity, R_d=0, size_file=ringfile)
WHEN(Para) AT(0,0,0) RELATIVE ref_p
ROTATED (0.0,0.0,0.0) RELATIVE ref_p
EXTEND
%{
    if (SCATTERED){
        parascatter=SCATTERED;
        pararef=p/pstore;
    }
%}

COMPONENT misalign_rot_h = Arm()
AT(dHx,dHy+PRM[ROW[porenumber-1]-1],dHz) RELATIVE a_1
ROTATED (dHrx/3600.0,dHry/3600.0,dHrz/3600.0) RELATIVE a_1

COMPONENT ref_h = Arm()
AT( 0,-PRM[ROW[porenumber-1]-1],0) RELATIVE misalign_rot_h

COMPONENT mm_h_1 = MM_h(
        pore_th=0, ring_nr=ROW[porenumber-1], Z0=FL, pore_width=pore_width, pore_height=pore_height, mirror_reflec=reflectivity, R_d=0, size_file=ringfile)
WHEN (Hyper) AT(0,0,0) RELATIVE ref_h
ROTATED (0.0,0.0,0.0) RELATIVE ref_h
EXTEND
%{
    if (SCATTERED){
        hyperscatter=SCATTERED;
        hyperref=p/(pstore*pararef);
    }
%}

COMPONENT detector_post_optics = PSD_monitor(restore_xray=1,xwidth=3, yheight=3, nx=101, 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_pore_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_pore_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