[neutron-mc] McStas Polarisation details

[Dalgliesh, RM (Robert)]

Hi Klaus, Peter, Kristian,
Ok, I've done more testing with a series simple components that report
the polarisation direction of the neutron and nothing else. After an
hour or so of head scratching thankfully Klaus is right. Using simple
components my main concern regarding what happens with absolute
positioning appears not to be a problem. The polarisation matrix is only
changed locally within each component and so the test instrument below
does exactly as it should.
i.e. if you rotate the first 2 pol_test components they report the
polarisation direction along the y direction but the 3rd if left
unrotated does not.
It serves me right for diving in at the deep end. I will have to
carefully check through my polariser component to work out exactly where
I was going wrong with my polarisation definitions. Hopefully I will
have a fully tested version soon. 

Thanks for the help and suggestions.

Rob

---------------------------------------------------
/***********************************************************************
********
*
* Component: pol_set
*
* %I
* Written by: Robert Dalgliesh
* Date: August 2005
* Version: $Revision: 0.5 $
* Origin: ISIS
*
* Component to define the polarisation of a neutorn beam
*
* %D
* The routine defines the polarisation of a neutron beam
* to be 50% +/- along either the x,y or z axes.
*
* Example: set_pol(sx1=1,sy1=0,sz1=0)
*
* %P
* sx1:     (1) Direction beam is to be polarised along x axis
* sy1:     (1) Direction beam is to be polarised along y axis
* sz1:     (1) Direction beam is to be polarised along z axis
*
* %E
************************************************************************
*******/
DEFINE COMPONENT pol_set
DEFINITION PARAMETERS ()
SETTING PARAMETERS (sx1=1, sy1=0, sz1=0)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
POLARISATION PARAMETERS (sx,sy,sz)

/* define random polarisation direction this can probably be done more
quickly*/

TRACE
%{
  double norm,sumpol;
  
  sumpol=sx1+sy1+sz1;
  if (sumpol > 1) 
  { fprintf(stderr,"pol_set: You have tried to define polarisation in
more than one direction\n",NAME_CURRENT_COMP); 
    exit(-1); }
  
  if(sx1 > 0){
    do{
      sx = randpm1();
      norm = sx*sx;
    } while (norm > 1 || norm == 0);
    norm=sqrt(norm);
    sx /= norm;
    sy = 0.0;
    sz = 0.0;
  }else if(sy1 > 0){
    do{
      sy = randpm1();
      norm = sy*sy;
    } while (norm > 1 || norm == 0);
    norm=sqrt(norm);
    sy /= norm;
    sx = 0.0;
    sz = 0.0;
  }else{
    do{
      sz = randpm1();
      norm = sz*sz;
    } while (norm > 1 || norm == 0);
    norm=sqrt(norm);
    sz /= norm;
    sx = 0.0;
    sy = 0.0;
  }
  printf("pol_set: sx %f, sy %f, sz %f \n",sx,sy,sz);
%}

END
-------------------------------------------------------
/***********************************************************************
********
*
* Component: pol_test
*
* %I
* Written by: Robert Dalgliesh
* Date: August 2005
* Version: $Revision: 0.1 $
* Origin: ISIS
*
* Do nothing
*
* %D
* do nothing but print polarisation direction
*
* Example: set_test
*
* %P
*
* %E
************************************************************************
*******/
DEFINE COMPONENT pol_test
DEFINITION PARAMETERS ()
SETTING PARAMETERS ()
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
POLARISATION PARAMETERS (sx,sy,sz)

TRACE
%{
 if(abs(sy) > 1e-6){
   printf("Pol along sy\n");
 }
 printf("pol_test: sx %f, sy %f, sz %f \n",sx,sy,sz);
%}

END
------------------------------------------------------
DEFINE INSTRUMENT test(az1=0.0,az2=0.0,az3=0.0)
DECLARE
%{
%}

INITIALIZE
%{
%}


TRACE

COMPONENT Origin = Arm()
  AT (0,0,0) ABSOLUTE

COMPONENT source1 = Source_flat_lambda(
    radius = 0.05, dist = 18.0, xw = 0.04, yh = 0.002,
    lambda_0 = 8.0, d_lambda = 7.5)
  AT (0.0, 0.0, 0.00001) RELATIVE Origin

COMPONENT pol_set1 = pol_set(
    sx1 = 1, sy1 = 0, sz1 = 0)
  AT (0.0, 0.0, 0.01) RELATIVE source1

COMPONENT slit1 = Slit(
    xmin = -0.02, xmax = 0.02, ymin = -0.001, ymax = 0.001)
  AT (0.0, 0.0, 18.0) RELATIVE source1

COMPONENT ptest1 = pol_test(
    )
  AT (0.0, 0.0, 1.0) RELATIVE slit1
  ROTATED (0.0, 0.0, az1) RELATIVE slit1

COMPONENT ptest2 = pol_test(
    )
  AT (0.0, 0.0, 0.1) RELATIVE ptest1
  ROTATED (0.0, 0.0, az2) RELATIVE ptest1

COMPONENT ptest3 = pol_test(
    )
  AT (0.0, 0.0, 1.2) RELATIVE slit1
  ROTATED (0.0, 0.0, az3) RELATIVE slit1

FINALLY
%{
%}
END