Add TOFtrain version of MCPL instr

Author: willend

mcstas-comps/examples/Prototypes/ODIN_MCPL_TOF_train4/DiskChopper.comp

@@ -0,0 +1,230 @@
+/*******************************************************************************
+*
+* McStas, neutron ray-tracing package
+*         Copyright (C) 1997-2008, All rights reserved
+*         Risoe National Laboratory, Roskilde, Denmark
+*         Institut Laue Langevin, Grenoble, France
+*
+* Component: DiskChopper
+*
+* %I
+* Written by: Peter Willendrup
+* Date: March 9 2006
+* Origin: Risoe
+* Based on Chopper (Philipp Bernhardt), Jitter and beamstop from work by
+* Kaspar Hewitt Klenoe (jan 2006), adjustments by Rob Bewey (march 2006)
+*
+* %D
+* Models a disc chopper with nslit identical slits, which are symmetrically distributed
+* on the disc. At time t=0, the centre of the first slit opening will be situated at the
+* vertical axis when phase=0, assuming the chopper centre of rotation is placed <b>BELOW</b> the beam axis.
+* If you want to place the chopper <b>ABOVE</b> the beam axis, please use a 180 degree rotation around Z 
+* (otherwise unexpected beam splitting can occur in combination with the isfirst=1 setting, see
+*  <a href="https://github.com/mccode-dev/McCode/issues/650">related bug on GitHub</a>)
+*
+* For more complicated gemometries, see component manual example of DiskChopper GROUPing.
+*
+* If the chopper is the 1st chopper of a continuous source instrument, you should use the "isfirst" parameter.
+* This parameter SETS the neutron time to match the passage of the chooper slit(s), taking into account the
+* chopper timing and phasing (thus conserving your simulated statistics).
+*
+* The isfirst parameter is ONLY relevant for use in continuous source settings.
+*
+* Example: DiskChopper(radius=0.2, theta_0=10, nu=41.7, nslit=3, delay=0, isfirst=1) First chopper
+*          DiskChopper(radius=0.2, theta_0=10, nu=41.7, nslit=3, delay=0, isfirst=0)
+*
+* NOTA BENE wrt. GROUPing and isfirst: 
+* When setting up a GROUP of DiskChoppers for a steady-state / reactor source, you will need
+* to set up
+* 1) An initial chopper with isfirst=1, NOT part of the GROUP - and using a "big" chopper opening 
+*    that spans the full angular extent of the openings of the subsequent GROUP
+* 2) Add your DiskChopper GROUP setting isfirst=0
+*
+* %P
+* INPUT PARAMETERS:
+*
+* theta_0: [deg]  Angular width of the slits.
+* yheight: [m]    Slit height (if = 0, equal to radius). Auto centering of beam at half height.
+* radius: [m]     Radius of the disc
+* nu: [Hz]        Frequency of the Chopper, omega=2*PI*nu (algebraic sign defines the direction of rotation)
+* nslit: [1]      Number of slits, regularly arranged around the disk
+*
+* Optional parameters:
+* isfirst: [0/1]  Set it to 1 for the first chopper position in a cw source (it then spreads the neutron time distribution)
+* n_pulse: [1]    Number of pulses (Only if isfirst)
+* jitter: [s]     Jitter in the time phase
+* abs_out: [0/1]  Absorb neutrons hitting outside of chopper radius?
+* delay: [s]      Time 'delay'
+* phase: [deg]    Angular 'delay' (overrides delay)
+* xwidth: [m]     Horizontal slit width opening at beam center
+* verbose: [1]    Set to 1 to display Disk chopper configuration
+*
+* %E
+*******************************************************************************/
+
+DEFINE COMPONENT DiskChopper
+
+
+
+SETTING PARAMETERS (theta_0=0, radius=0.5, yheight, nu, nslit=3, jitter=0, delay=0, isfirst=0, n_pulse=1, abs_out=1, phase=0, xwidth=0, verbose=0)
+
+
+/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
+
+DECLARE
+%{
+  double Tg;
+  double To;
+  double delta_y;
+  double height;
+  double omega;
+%}
+
+INITIALIZE
+%{
+  /* If slit height 'unset', assume full opening */
+  if (yheight == 0) {
+    height = radius;
+  } else {
+    height = yheight;
+  }
+  delta_y = radius - height / 2; /* radius at beam center */
+  omega = 2.0 * PI * nu;         /* rad/s */
+  if (xwidth && !theta_0 && radius)
+    theta_0 = 2 * RAD2DEG * asin (xwidth / 2 / delta_y);
+
+  if (nslit <= 0 || theta_0 <= 0 || radius <= 0) {
+    fprintf (stderr, "DiskChopper: %s: nslit, theta_0 and radius must be > 0\n", NAME_CURRENT_COMP);
+    exit (-1);
+  }
+  if (nslit * theta_0 >= 360) {
+    fprintf (stderr, "DiskChopper: %s: nslit * theta_0 exceeds 2PI\n", NAME_CURRENT_COMP);
+    exit (-1);
+  }
+  if (yheight && yheight > radius) {
+    fprintf (stderr, "DiskChopper: %s: yheight must be < radius\n", NAME_CURRENT_COMP);
+    exit (-1);
+  }
+  if (isfirst && n_pulse <= 0) {
+    fprintf (stderr, "DiskChopper: %s: wrong First chopper pulse number (n_pulse=%g)\n", NAME_CURRENT_COMP, n_pulse);
+    exit (-1);
+  }
+  if (!omega) {
+    fprintf (stderr, "DiskChopper: %s WARNING: chopper frequency is 0!\n", NAME_CURRENT_COMP);
+    omega = 1e-15; /* We should actually use machine epsilon here... */
+  }
+  if (!abs_out) {
+    fprintf (stderr, "DiskChopper: %s WARNING: chopper will NOT absorb neutrons outside radius %g [m]\n", NAME_CURRENT_COMP, radius);
+  }
+
+  theta_0 *= DEG2RAD;
+
+  /* Calulate delay from phase and vice versa */
+  if (phase) {
+    if (delay) {
+      fprintf (stderr, "DiskChopper: %s WARNING: delay AND phase specified. Using phase setting\n", NAME_CURRENT_COMP);
+    }
+    phase *= DEG2RAD;
+    /* 'Delay' should always be a delay, taking rotation direction into account: */
+    delay = phase / fabs (omega);
+  } else {
+    phase = delay * omega; /* rad */
+  }
+
+  /* Time from opening of slit to next opening of slit */
+  Tg = 2.0 * PI / fabs (omega) / nslit;
+
+  /* How long can neutrons pass the Chopper at a single point */
+  To = theta_0 / fabs (omega);
+
+  if (!xwidth)
+    xwidth = 2 * delta_y * sin (theta_0 / 2);
+
+  if (verbose && nu) {
+    printf ("DiskChopper: %s: frequency=%g [Hz] %g [rpm], time frame=%g [s] phase=%g [deg]\n", NAME_CURRENT_COMP, nu, nu * 60, Tg, phase * RAD2DEG);
+    printf ("             %g slits, angle=%g [deg] height=%g [m], width=%g [m] at radius=%g [m]\n", nslit, theta_0 * RAD2DEG, height, xwidth, delta_y);
+  }
+%}
+
+TRACE
+%{
+  double toff;
+  double yprime;
+  PROP_Z0;
+  yprime = y + delta_y;
+
+  /* Is neutron outside the vertical slit range and should we absorb? */
+  if (abs_out && (x * x + yprime * yprime) > radius * radius) {
+    ABSORB;
+  }
+  /* Does neutron hit inner solid part of chopper in case of yheight!=radius? */
+  if ((x * x + yprime * yprime) < (radius - height) * (radius - height)) {
+    ABSORB;
+  }
+
+  if (isfirst) {
+    /* all events are put in the transmitted time frame */
+    t = atan2 (x, yprime) / omega + To * randpm1 () / 2.0 + delay + (jitter ? jitter * randnorm () : 0) + (n_pulse > 1 ? floor (n_pulse * rand01 ()) * Tg : 0);
+    /* correction: chopper slits transmission opening/full disk */
+    p *= nslit * theta_0 / 2.0 / PI;
+  } else {
+	  
+	// Check whether each t_offset carried by the ray make it through
+	double weight_update = p/_particle->p_last_time_manipulation;
+	_particle->p_last_time_manipulation = 0;
+	
+    int train_index;
+	int one_did_hit = 0; 
+	double this_train_t;
+	int all_dead = 1;
+
+    for (train_index=0; train_index<adaptive_N; train_index++) {
+		
+		if (_particle->p_trains[train_index] == 0) continue;
+		all_dead = 0;
+		
+		this_train_t = t + _particle->t_offset[train_index];
+        toff = fabs (this_train_t - atan2 (x, yprime) / omega - delay - (jitter ? jitter * randnorm () : 0));
+
+        /* does neutron hit outside slit? */
+        if (fmod (toff + To / 2.0, Tg) > To)
+          _particle->p_trains[train_index] = 0; // T_ABSORB
+		else {	
+			// T_TRANSMIT
+		   one_did_hit = 1;
+		   _particle->p_trains[train_index] *= weight_update;
+		   _particle->p_last_time_manipulation += _particle->p_trains[train_index];
+	   }
+		
+	}
+	if (!one_did_hit || all_dead) ABSORB;
+	
+	p = _particle->p_last_time_manipulation;
+	
+  }
+  SCATTER;
+%}
+
+MCDISPLAY
+%{
+
+  int j;
+  /* Arrays for storing geometry of slit/beamstop */
+
+  circle ("xy", 0, -delta_y, 0, radius);
+
+  /* Drawing the slit(s) */
+  for (j = 0; j < nslit; j++) {
+    /* Angular start/end of slit */
+    double tmin = j * (2.0 * PI / nslit) - theta_0 / 2.0 + phase;
+    double tmax = tmin + theta_0;
+    /* Draw lines for each slit. */
+
+    line (radius * sin (tmin), radius * cos (tmin) - delta_y, 0, (radius - height) * sin (tmin), (radius - height) * cos (tmin) - delta_y, 0);
+    line ((radius - height) * sin (tmin), (radius - height) * cos (tmin) - delta_y, 0, (radius - height) * sin (tmax), (radius - height) * cos (tmax) - delta_y,
+          0);
+    line ((radius - height) * sin (tmax), (radius - height) * cos (tmax) - delta_y, 0, radius * sin (tmax), radius * cos (tmax) - delta_y, 0);
+  }
+%}
+
+END