McStas: forgotten attachment

[Pay H. Schulze Horn]

/* This Monitor give the horizontal divergence in korrelation to the energy in a 2D McStas styled array out.
If you want the vertikal divergence, turn it like you will.

THERE IS NO SECTION FOR MCDISPLAY!

* xmin:     Lower x bound of detector opening (m)
* xmax:     Upper x bound of detector opening (m)
* ymin:     Lower y bound of detector opening (m)
* ymax:     Upper y bound of detector opening (m)
* nenergy:  Number of counting-rate-channels for the intervall [e_min;e_max] of the energy.
* ndiv:		Number of counting-rate-channels for the intervall [-maxdiv;maxdiv] of the horizontal divergence.
* maxdiv:	Maximum divergence to detect.
* E_middle: The middle of the energy intervall [E_middle-E_delta;E_middle+E_delta] to detect.
* E_delta: 	Gives the range of the energy intervall [E_middle-E_delta;E_middle+E_delta] to detect.

This Monitor is component for McStas and written by Pay Schulze Horn at 
Institut Laue Langevin in August 1999.
e-mail: horn at ill.fr until 26.10.99
		  stu32232 at mail.uni-kiel.de for the next 2 years

McStas is written by Kim Lefmann and
Kristian Nielsen at Risoe National Laboratory Danmark.

The source-code is based on the component-source-code Divergence_monitor, which
is written by: KL,  Nov. 11, 1998 

*/

DEFINE COMPONENT Divenergymon
DEFINITION PARAMETERS (xmin, xmax, ymin, ymax, nenergy, ndiv, maxdiv,E_middle,E_delta,filename)
SETTING PARAMETERS ()
OUTPUT PARAMETERS (Div_N, Div_p, Div_p2) 
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
  %{
    
	 /* the 2D-arrays in which the data will be stored */
	 int Div_N[ndiv][nenergy];
    double Div_p[ndiv][nenergy];
    double Div_p2[ndiv][nenergy];
  %}
INITIALIZE
  %{
	 int i,j;
    for (i=0; i<ndiv; i++)
     for (j=0; j<nenergy; j++) 
     {
      Div_N[i][j] = 0;
      Div_p[i][j] = 0;
      Div_p2[i][j] = 0;
     }
  %}
TRACE
  %{
    int i,j;
    double neut_energy, neut_div;
    /* Here are getting all neutrons lost, which won't fly to the z=0 plane. */
	 PROP_Z0;
    
	 /* Only neutrons that will fly through the monitor area will be rocognized. */
	 if (x>xmin && x<xmax && y>ymin && y<ymax)
  	 {
		neut_energy =(vx*vx+vy*vy+vz*vz)*VS2E;
		neut_div = RAD2DEG*atan2(vx,vz);
      
		/* making sure that the neutron is in the described intervalls */
		if (neut_energy < E_middle+E_delta && neut_energy > E_middle-E_delta && 
          neut_div < maxdiv && neut_div > -maxdiv)
      {
        /* counting at what space the datas will be stored */
		  i = floor(((neut_energy-E_middle+E_delta)/(2*E_delta))*nenergy);
        j = floor((neut_div + maxdiv)*ndiv/(2.0*maxdiv));
        Div_N[i][j]++;
        Div_p[i][j] += p;
        Div_p2[i][j] += p*p;
      }
    }
  %}
FINALLY
  %{
    DETECTOR_OUT_2D(
        "horizontal divergence related to energy monitor",
        "Energy [meV]",
        "horizontal(y) divergence [deg]",
        E_middle-E_delta, E_middle+E_delta, -maxdiv, maxdiv,
        ndiv, nenergy,
        &Div_N[0][0],&Div_p[0][0],&Div_p2[0][0],
        filename);
  %}

END