/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2006, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: SiC
*
*%I
* Written by: <a href="mailto:S.RYCROFT@IRI.TUDELFT.NL">S. Rycroft</a>
* Date: Jan 2000
* Origin: IRI.
* Release: McStas 1.9.1
*
* SiC layer sample
*
* %D
* SiC layer sample
* 
* Example: SiC()
*
* %P
* INPUT PARAMETERS:
*
* xlenght: (m)    length of SiC plate 
* yheight: (m)    height of SiC plate
*
* %E
*******************************************************************************/



DEFINE COMPONENT SiC
DEFINITION PARAMETERS (xlength, yheight)
SETTING PARAMETERS ()
OUTPUT PARAMETERS(ScatGrad)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

DECLARE
%{
  double ScatGrad[85];
%}


INITIALIZE
%{
  int i;
  
  for (i = 0; i < 85; ++i)
    ScatGrad[i] = 0;
  ScatGrad[1]=7.3106e-6;
  ScatGrad[81]=2.073e-6-7.3106e-6;
%}


TRACE
%{
  double dt, q;

/* Variables added for Rayleigh Appproximation. */

  double realZ,imagZ;
  double csZ,snZ,Qc,R0;
  int i,Z,q_count;

  Qc=0.010208;
  R0=1.0;

  /* First check if neutron has the right direction. */

  if(vz != 0.0 && (dt = -z/vz) >= 0)
  {
    double old_x = x, old_y = y;

    x += vx*dt;
    y += vy*dt;
    /* Now check if neutron intersects mirror. */
    if(x >= 0 && x <= xlength && y >= 0 && y <= yheight)
    {
      z = 0;
      t += dt;
      q = fabs(2*vz*V2Q);
      vz = -vz;
      /* Reflectivity (see component Guide).
   	     Changed to calculate from real sample. */

      if(q > Qc)
      {
        realZ = 0;
        imagZ = 0;
        for (Z = 0; Z < 85; Z++)
        {
          csZ = cos(-1.*q*Z);
          snZ = sin(-1.*q*Z);
          realZ = realZ + ScatGrad[Z] * csZ;
          imagZ = imagZ + ScatGrad[Z] * snZ;
	      }
        p *=16.*PI*PI*((realZ*realZ)+(imagZ*imagZ))/(q*q*q*q);
      }
      p *= R0;
      SCATTER;
    }
    else
    {
      x = old_x;
      y = old_y;
    }
  }
%}

END








