/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Instrument: Prisma2
*
* %Identification
* Written by: Kristian Nielsen and Mark Hagen
* Date: August 1998.
* Origin: ISIS/Risoe
* Release: McStas 1.1
* Version: $Revision: 1.6 $
* %INSTRUMENT_SITE: ISIS
*
* Simple simulation of PRISMA2 with RITA-style analyser backend.
*
* %Description
* Demonstrates how the standard components from the component library
* may be easily modified for special purposes; in this case to have
* the individual analyser blades paint a "color" on the neutrons to
* differentiate them in the detector.
*
* Output is in the file "prisma2.tof". The format is ASCII; each
* line consists of the time-of-flight in microseconds followed by seven
* intensities of neutrons from each individual analyser blade.
*
* Example: mcrun prisma2.instr TT=-30
*
* %Parameters
* INPUT PARAMETERS:
* TT:  (deg) Take-off angle at the sample position, aka A4
* PHA: (deg) Analyzer group rotation angle, aka A5
* PHA1:(deg) Analyzer 1 tilt angle
* PHA2:(deg) Analyzer 2 tilt angle
* PHA3:(deg) Analyzer 3 tilt angle
* PHA4:(deg) Analyzer 4 tilt angle
* PHA5:(deg) Analyzer 5 tilt angle
* PHA6:(deg) Analyzer 6 tilt angle
* PHA7:(deg) Analyzer 7 tilt angle
* TTA: (deg) Take-off angle at the analyzer position, aka A6
*
* %Link
* The McStas User manual
* <a href="http://www.isis.rl.ac.uk/excitations/prisma/">PRISMA</a>
*
* %End
*******************************************************************************/

/* Modified from Monochromator.comp to paint a "color" on the neutron
   if it is scattered. */
DEFINE COMPONENT Monochromator_color
DEFINITION PARAMETERS (zmin, zmax, ymin, ymax, mosaich, mosaicv, r0, Q, color)
SETTING PARAMETERS ()
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
  %{
#define DIV_CUTOFF 2            /* ~ 10^-5 cutoff. */
  %}

TRACE
  %{
    double dphi,tmp1,tmp2,tmp3,vratio,phi,theta0,theta,v,cs,sn;
    double old_x = x, old_y = y, old_z = z, old_t = t;
    double dt;

    if(vx != 0.0 && (dt = -x/vx) >= 0.0)
    {
      y += vy*dt; z += vz*dt; t += dt; x = 0.0;

    if (z>zmin && z<zmax && y>ymin && y<ymax)
    {
      /* First: scattering in plane */

      theta0 = atan2(vx,vz);           /* neutron angle to slab */
      v = sqrt(vx*vx+vy*vy+vz*vz);
      theta = asin(Q2V*Q/(2.0*v));               /* Bragg's law */
      if(theta0 < 0)
        theta = -theta;
      tmp3 = (theta-theta0)/(MIN2RAD*mosaich);
      if(tmp3 > DIV_CUTOFF)
      {
        x = old_x; y = old_y; z = old_z; t = old_t;
      }
      else
      {
        p *= r0*exp(-tmp3*tmp3*4*log(2)); /* Use mosaics */
        tmp1 = 2*theta;
        cs = cos(tmp1);
        sn = sin(tmp1);
        tmp2 = cs*vx - sn*vz;
        vy = vy;
        vz = cs*vz + sn*vx;
        vx = tmp2;

        /* Second: scatering out of plane.
           Approximation is that Debye-Scherrer cone is a plane */

        phi = atan2(vy,vz);                            /* out-of plane angle */
        dphi = (MIN2RAD*mosaicv)/(2*sqrt(2*log(2)))*randnorm();  /* MC choice: */
        /* Vertical angle of the crystallite */
        vy = vz*tan(phi+2*dphi*sin(theta));
        vratio = v/sqrt(vx*vx+vy*vy+vz*vz);
        vz = vz*vratio;
        vy = vy*vratio;                             /* Renormalize v */
        vx = vx*vratio;
        neu_color = color;
      }
    }
    else
    {
      x = old_x; y = old_y; z = old_z; t = old_t;
    }
    }
  %}
MCDISPLAY
  %{
    magnify("zy");
    multiline(5, 0.0, (double)ymin, (double)zmin,
                 0.0, (double)ymax, (double)zmin,
                 0.0, (double)ymax, (double)zmax,
                 0.0, (double)ymin, (double)zmax,
                 0.0, (double)ymin, (double)zmin);
  %}
END


/* Modified from TOF_monitor.comp to bin neutrons according to their
   "color". */
DEFINE COMPONENT TOF_monitor_color
DEFINITION PARAMETERS (xmin, xmax, ymin, ymax, nchan, dt, filename, maxcolor)
SETTING PARAMETERS ()
OUTPUT PARAMETERS (TOF_N, TOF_p, TOF_p2)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
  %{
#ifdef MCSTAS_VERSION
    double TOF_N[maxcolor+1][nchan];
#else
    int    TOF_N[maxcolor+1][nchan];
#endif
    double TOF_p[maxcolor+1][nchan];
    double TOF_p2[maxcolor+1][nchan];
  %}
INITIALIZE
  %{
    int i,c;

    for (i=0; i<nchan; i++)
      for (c=0; c<=maxcolor; c++)
      {
        TOF_N[c][i] = 0;
        TOF_p[c][i] = 0;
        TOF_p2[c][i] = 0;
      }
  %}
TRACE
  %{
    int i;

    PROP_Z0;
    if (x>xmin && x<xmax && y>ymin && y<ymax)
    {
      i = floor(1E6*t/dt);              /* Bin number */
      if(i >= nchan) i = nchan;
      if(i < 0)
      {
        printf("FATAL ERROR: negative time-of-flight.\n");
        exit(1);
      }
      if(neu_color < 0 || neu_color > maxcolor)
      {
        printf("FATAL ERROR: wrong color neutron.\n");
        exit(1);
      }
      TOF_N[neu_color][i]++;
      TOF_p[neu_color][i] += p;
      TOF_p2[neu_color][i] += p*p;
    }
  %}
FINALLY
  %{
    printf("writing TOF color\n");

    DETECTOR_OUT_2D(
        "Time-of-flight monitor",
        "Neutron color",
        "Time-of-flight [us]",
        0, maxcolor+1, 0, nchan*10,
        maxcolor+1, nchan,
        &TOF_N[0][0],&TOF_p[0][0],&TOF_p2[0][0],
        filename);

  %}

MCDISPLAY
  %{
    magnify("xy");
    multiline(5, (double)xmin, (double)ymin, 0.0,
                 (double)xmax, (double)ymin, 0.0,
                 (double)xmax, (double)ymax, 0.0,
                 (double)xmin, (double)ymax, 0.0,
                 (double)xmin, (double)ymin, 0.0);
  %}
END


DEFINE INSTRUMENT prisma2(TT=-30,PHA=22,PHA1=-3,PHA2=-2,PHA3=-1,
                          PHA4=0,PHA5=1,PHA6=2,PHA7=3,TTA=4)

DECLARE
%{
int neu_color;                  /* "Color" of current neutron */

/* 30' mosaicity used on analysator */
double prisma_ana_mosaic = 30;
/* Q vector for bragg scattering with monochromator and analysator */
double prisma_ana_q = 1.87325;
double prisma_ana_r0 = 0.6;
double focus_x,focus_z;

double apos1, apos2, apos3, apos4, apos5, apos6, apos7;
%}

INITIALIZE
%{
  focus_x = 0.52 * sin(TT*DEG2RAD);
  focus_z = 0.52 * cos(TT*DEG2RAD);
  /* Rita-style analyser. */
  {
    double l = 0.0125;
    apos1 = -3*l;
    apos2 = -2*l;
    apos3 = -1*l;
    apos4 =  0*l;
    apos5 =  1*l;
    apos6 =  2*l;
    apos7 =  3*l;
  }
%}

TRACE

COMPONENT mod = Moderator(
                          radius = 0.0707,
                          dist = 9.035,
                          xw = 0.021,
                          yh = 0.021,
                          E0 = 10, E1 = 15,
                          Ec = 9.0, t0 = 37.15, gam = 39.1)
     AT (0,0,0) ABSOLUTE

/* Use a slit to get the effect of a rectangular source. */
COMPONENT modslit = Slit(xmin = -0.05, xmax = 0.05,
                         ymin = -0.05, ymax = 0.05)
     AT(0,0,0.000001) RELATIVE mod

COMPONENT tof_test = TOF_monitor(xmin = -0.05, xmax = 0.05,
                                 ymin = -0.05, ymax = 0.05,
                                 nchan = 500, dt = 1,
                                 filename = "prisma2.mon")
     AT (0,0,0.005) RELATIVE mod

COMPONENT mon1 = Monitor(xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
     AT(0,0,0.01) RELATIVE mod ROTATED (0,0,0) RELATIVE mod

COMPONENT slit1 = Slit(xmin = -0.05, xmax = 0.05,
                         ymin = -0.05, ymax = 0.05)
     AT(0,0,1.7) RELATIVE mod

COMPONENT slit2 = Slit(xmin = -0.02, xmax = 0.02,
                         ymin = -0.03, ymax = 0.03)
     AT(0,0,7) RELATIVE slit1

COMPONENT mon2 = Monitor(xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
     AT(0,0,9) RELATIVE mod

COMPONENT sample = V_sample(
        radius_i = 0.00001, radius_o = 0.01,
        h = 0.02,
        focus_r = 0.03,
        pack = 1,
        target_x = focus_x, target_y = 0, target_z = focus_z)
  AT (0, 0, 9.035) RELATIVE mod

COMPONENT a2 = Arm() AT (0,0,0) RELATIVE sample ROTATED (0,TT,0) RELATIVE sample

COMPONENT mon3 = Monitor(xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
     AT(0,0,0.39) RELATIVE a2

COMPONENT coll2 = Soller(xmin = -0.015, xmax = 0.015,
                         ymin = -0.025, ymax = 0.025,
                         len = 0.12, divergence = 120)
     AT(0,0,0.40) RELATIVE a2

COMPONENT mon4 = Monitor(xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
     AT(0,0,0.521) RELATIVE a2

COMPONENT rita_ana = Arm()
     AT(0, 0, 0.58) relative a2 ROTATED (0, PHA, 0) RELATIVE a2

COMPONENT ana1 = Monochromator_color(
        ymin=-0.0375,ymax=0.0375,zmin=-0.006,zmax=0.006,
        mosaich=prisma_ana_mosaic,mosaicv=prisma_ana_mosaic,
        r0=prisma_ana_r0, Q=prisma_ana_q, color = 0)
  AT (0, 0, apos1) RELATIVE rita_ana
  ROTATED (0, PHA1, 0) RELATIVE rita_ana

COMPONENT ana2 = Monochromator_color(
        ymin=-0.0375,ymax=0.0375,zmin=-0.006,zmax=0.006,
        mosaich=prisma_ana_mosaic,mosaicv=prisma_ana_mosaic,
        r0=prisma_ana_r0, Q=prisma_ana_q, color = 1)
  AT (0, 0, apos2) RELATIVE rita_ana
  ROTATED (0, PHA2, 0) RELATIVE rita_ana

COMPONENT ana3 = Monochromator_color(
        ymin=-0.0375,ymax=0.0375,zmin=-0.006,zmax=0.006,
        mosaich=prisma_ana_mosaic,mosaicv=prisma_ana_mosaic,
        r0=prisma_ana_r0, Q=prisma_ana_q, color = 2)
  AT (0, 0, apos3) RELATIVE rita_ana
  ROTATED (0, PHA3, 0) RELATIVE rita_ana

COMPONENT ana4 = Monochromator_color(
        ymin=-0.0375,ymax=0.0375,zmin=-0.006,zmax=0.006,
        mosaich=prisma_ana_mosaic,mosaicv=prisma_ana_mosaic,
        r0=prisma_ana_r0, Q=prisma_ana_q, color = 3)
  AT (0, 0, apos4) RELATIVE rita_ana
  ROTATED (0, PHA4, 0) RELATIVE rita_ana

COMPONENT ana5 = Monochromator_color(
        ymin=-0.0375,ymax=0.0375,zmin=-0.006,zmax=0.006,
        mosaich=prisma_ana_mosaic,mosaicv=prisma_ana_mosaic,
        r0=prisma_ana_r0, Q=prisma_ana_q, color = 4)
  AT (0, 0, apos5) RELATIVE rita_ana
  ROTATED (0, PHA5, 0) RELATIVE rita_ana

COMPONENT ana6 = Monochromator_color(
        ymin=-0.0375,ymax=0.0375,zmin=-0.006,zmax=0.006,
        mosaich=prisma_ana_mosaic,mosaicv=prisma_ana_mosaic,
        r0=prisma_ana_r0, Q=prisma_ana_q, color = 5)
  AT (0, 0, apos6) RELATIVE rita_ana
  ROTATED (0, PHA6, 0) RELATIVE rita_ana

COMPONENT ana7 = Monochromator_color(
        ymin=-0.0375,ymax=0.0375,zmin=-0.006,zmax=0.006,
        mosaich=prisma_ana_mosaic,mosaicv=prisma_ana_mosaic,
        r0=prisma_ana_r0, Q=prisma_ana_q, color = 6)
  AT (0, 0, apos7) RELATIVE rita_ana
  ROTATED (0, PHA7, 0) RELATIVE rita_ana


COMPONENT a3 = Arm()
     AT (0,0,0) relative rita_ana ROTATED (0,TTA,0) RELATIVE a2

COMPONENT mon5 = Monitor(xmin = -0.05, xmax = 0.05, ymin = -0.05, ymax = 0.05)
     AT(0,0,0.06) RELATIVE a3

COMPONENT mon6 = Monitor(xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
     AT(0,0,0.161) RELATIVE a3

COMPONENT psd = PSD_monitor(xmin = -0.05, xmax = 0.05,
                            ymin = -0.05, ymax = 0.05,
                            nx = 100, ny = 100,
                            filename = "prisma2.psd")
     AT(0,0,0.20) RELATIVE a3

COMPONENT detector = TOF_monitor_color(xmin = -0.05, xmax = 0.05,
                                 ymin = -0.05, ymax = 0.05,
                                 nchan = 10000, dt = 10, maxcolor = 6,
                                 filename = "prisma2.tof")
     AT (0,0,0.20) RELATIVE a3

COMPONENT mon9 = Monitor(xmin = -0.1, xmax = 0.1, ymin = -0.1, ymax = 0.1)
     AT(0,0,0.01) RELATIVE detector

END