/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: my_component
*
* %IDENTIFICATION
* Written by: Erik B Knudsen
* Date: Sep. 2015
* Origin: DTU Physics
* Modified by: 
*
* %DESCRIPTION
* A magnet with a spin-flip foil inside
* 
* This component models a magnet with a constant magnetic field and a slanted 
* magnetized foil inside acting as a spin-flipper. This arrangement may be used to
* target the geometry of a Spin-Echo SANS instrument (Rekveldt, NIMB, 1997 || Rekveldt, Rev Sci Instr. 2005).
* In its most basic form a SE-SANS instrument nrelies  on inclined field regions of
* constant magnetic field, which can be difficult to realize. Instead it is possible to emulate
* such regions using a magnetized foil.
*
* In this component th emagntzied foil is modelled as a canted (by the angle phi) mathematical plane inside a region
* of constant magnetic field of dimensions (xwidth,yheight,zdepth). Furthermore, exponetially decaying stray fields from
* the constant field region may be specified by giving parameters (Bxwidth,Byheight,Bzdepth) > (xwidth,yheight,zdepth).
*
*
* %PARAMETERS
* xwidth: [m]    width of the magnet 
* yheight: [m]   height of the magnet 
* zdepth: [m]    thickness of the magnet 
* Bxwidth: [m]   width of the B-field of thev component 
* Byheight: [m]  height of the B-field of the component 
* Bzdepth: [m]   thickness of the B-field of the component 
* phi: [rd]      angle (from the xz-plane) of the foil-spinflipper inside the magnet 
* foilthick: [um] Thickness of the magnetized foil
* Bx: [T]        Parameter used for x composant of field.
* By: [T]        Parameter used for y composant of field.
* Bz: [T]        Parameter used for z composant of field.
* stray_field: [ ] Toggle for modelling of stray fields.
* verbose: [ ]   Toggle verbose mode
* foil_in: [ ]   Toggle for optional foil in the flipper.
*
* %LINKS
*
* %END
****************************************************************************/

/*
This comment will not be exported to mcdoc
*/

DEFINE COMPONENT Foil_flipper_magnet

SETTING PARAMETERS(stray_field=1,xwidth, yheight, zdepth, Bxwidth=-1, Byheight=-1, Bzdepth=-1, phi=0.0, foilthick=0.0, Bx,By,Bz,foil_in=1, verbose=0)

SHARE
%{
/* Declare structures and functions only once in each instrument. */
#ifndef POL_LIB_H
%include "pol-lib"
#endif

#pragma acc routine
void foil_spin_rot(double *sx, double *sy, double *sz, double *vx, double *vy, double *vz, double phi, double th, double b){
  /*the spin flip is actually a precession around the field vector - in this case this vector is
   * in the foil plane, perp. to X, i.e. Z rotated by phi.*/
  double ux,uy,uz,s,cx,cy,cz;
  double theta,v,lamb,teff;

  /*now use (0,0,1) as a vector to rotate (flip) around to mimic the use of correction coils.*/
  ux=uy=0;uz=1;
  //ux=0;uy=sin(phi);uz=cos(phi);
  if (b < 0) {
    uy = -uy;
    uz = -uz;
  }
  v=sqrt(*vx * *vx + *vy * *vy + *vz * *vz);
  lamb = 2*PI*K2V/v;
  teff = th/sin(phi);
  /* Calculate the flipping angle theta by c*t*Bs*lambda
   * where c is the neutron precession constant [T^-1*m^-2], t is the
   * foil thickness [um], Bs is the induced magnetic field strength
   * in the foil [T] and lambda is the neutron wavelength [m] */
  theta = 4.6368*1e14*teff*1e-6*1*lamb*1e-10;
  /*Rodigues formula for roatating a vector around a unit vector.
   * v_rot=v*cos(theta) + uxv * sin(theta) + (u dot v)u(1-cos(theta)*/

  cx=*sz*uy-*sy*uz;
  cy=*sx*uz-*sz*ux;
  cz=*sy*ux-*sx*uy;
  s=scalar_prod(ux,uy,uz,*sx,*sy,*sz);
  *sx=*sx*cos(theta) + cx*sin(theta) + s*ux*(1-cos(theta));
  *sy=*sy*cos(theta) + cy*sin(theta) + s*uy*(1-cos(theta));
  *sz=*sz*cos(theta) + cz*sin(theta) + s*uz*(1-cos(theta));
}

#pragma acc routine
void foil_spin_flip(double *sx, double *sy, double *sz, double lamb, double phi){
  /*the spin flip is actually a precession around the field vector - in this case this vector is
   * in the foil plane, perp. to X, i.e. Z rotated by phi.*/
  double ux,uy,uz,s;
  ux=0;uy=sin(phi);uz=cos(phi);
  /*Rodigues formula for roatating a vector around a unit vector.
   * v_rot=v*cos(theta) + uxv * sin(theta) + (u dot v)u(1-cos(theta)
   * theta=180 deg. => v_rot =-v + 2*(u dot v)u*/
  s=scalar_prod(ux,uy,uz,*sx,*sy,*sz);
  *sx=-*sx+2*s*ux;
  *sy=-*sy+2*s*uy;
  *sz=-*sz+2*s*uz;
}

int exp_dec_magnetic_field (double x, double y, double z, double t, double *Bx, double *By, double *Bz, void *data){
  double *prms = (double *)data;
  double b;
  /*parameters are in the order, Bx0,By0,Bz0,z0, where z0 is the point from which B is decaying*/
  if(z<prms[3]){
    b=1;
  }else{
    b=exp(- fabs(z-prms[3]));
  }
  *Bx=prms[0]*b;
  *By=prms[1]*b;
  *Bz=prms[2]*b;
}

%}

DECLARE
%{
  /*Declarations of variables used in the whole of the component*/
  void *magnet_prms;
  double foil_n[3];
  Rotation spin_flip;
  double inside_magnet[3];
  double instray_magnet[4];
  double outstray_magnet[4];
  int hit_foil;
%}

INITIALIZE
%{
  /*default is to have the foil horizontal*/
  Rotation foil_rot;
  Coords n,tmp;

  n.x=0;n.y=1;n.z=0;
  rot_set_rotation(foil_rot,phi,0,0);
  tmp=rot_apply(foil_rot,n);
  foil_n[0]=tmp.x;foil_n[1]=tmp.y; foil_n[2]=tmp.z;
  /*setup spin-flip rotation. This is equivalent to mirroring in the foil plane.*/
  if(verbose){
    printf("INFO: (%s) Foil normal: (%g %g %g)\n",NAME_CURRENT_COMP, foil_n[0],foil_n[1],foil_n[2]); 
  }
  /*copy-paste from simpleBfield.comp*/
  /*this should be a store_magnet type function in pol-lib i think*/
  Coords localG = rot_apply(ROT_A_CURRENT_COMP, coords_set(0,-GRAVITY,0)); 
  
  /*initialize magnet parameters*/
  instray_magnet[1]=By;
  instray_magnet[2]=Bz;
  instray_magnet[3]=-zdepth/2.0;
  
  outstray_magnet[0]=Bx;
  outstray_magnet[1]=By;
  outstray_magnet[2]=Bz;
  outstray_magnet[3]=zdepth/2.0;

  inside_magnet[0]=Bx;
  inside_magnet[1]=By;
  inside_magnet[2]=Bz;

  /*if the magnetic field dimensions are not given, we assume there are no stray fields.
    I.e. the magnetic field does not extend beyond the magnet gap*/
/*  if (Bxwidth!=-1 && Byheight!=-1 && Bzdepth!=-1){*/
/*    stray_field=1;*/
/*  }*/
  if (Bxwidth==-1) Bxwidth=xwidth;
  if (Byheight==-1) Byheight=yheight;
  if (Bzdepth==-1) Bzdepth=zdepth;
  if (verbose){
      printf("INFO (%s): xw,yh,zd=(%g %g %g), (Bxw,Byh,Bzd)=(%g %g %g)\n",NAME_CURRENT_COMP,xwidth,yheight,zdepth,Bxwidth,Byheight,Bzdepth);
      if (foil_in)  printf("INFO (%s): Foil is IN\n",NAME_CURRENT_COMP);
      else printf("INFO (%s): Foil is OUT\n", NAME_CURRENT_COMP);
  }
%}

TRACE
%{
  double t1,t2,t3,t4,dt;
  double sp,v;
  mcmagnet_field_info *old_magnet,*outer_magnet;
  double *dd;
  v=sqrt(vx*vx + vy*vy + vz*vz);
  /*check if we hit component at all*/
  if (plane_intersect(&t1,x,y,z,vx,vy,vz,0,0,1,0,0,-Bzdepth/2.0)==0){
      if(verbose) fprintf(stderr,"INFO: (%s) Missed the magnetic field plane.\n",NAME_CURRENT_COMP);
      ABSORB;
  }
  if ( (box_intersect(&t1, &t2, x, y, z, vx, vy, vz, Bxwidth, Byheight, Bzdepth))==0){
  /*propagate to the start of the component (which is -Bzdepth/2)*/
      if(verbose) fprintf(stderr,"INFO: (%s) Missed the magnetic field.\n",NAME_CURRENT_COMP);
      ABSORB;
  }
  if (t1<0){
      if(verbose) fprintf(stderr,"WARNING (%s): Neutron is already inside component on entry.\n",NAME_CURRENT_COMP);
  }else{
      PROP_DT(t1);
      t2=t2-t1;
  }
  /*There is an implicit assumption here that the neutron has entered at z<zd/2. This might not be true.*/
  if (stray_field){
    /*neutron is now at the start of our b-field description*/
    /*this would be the place to push the magnetic field description from the magnet stack*/
    /*inject our magnet description into the propagation and remember the old one*/
    mcmagnet_push(_particle,0,&(ROT_A_CURRENT_COMP),&(POS_A_CURRENT_COMP),0,instray_magnet);
  }
  /*do we hit the magnet opening assuming the magnet gap "tunnel" is rectangular*/
  if ( (box_intersect(&t3, &t4, x, y, z, vx, vy, vz,xwidth, yheight, zdepth))==0){
    if(verbose) fprintf(stdout,"INFO: (%s) Missed the magnet gap\n",NAME_CURRENT_COMP);
    ABSORB;
  }
  /* t3 is the time for intersecting the first part of the area under the magnet - propagate to that*/
  if (t3>1e-9){
    PROP_DT(t3);
    t4=t4-t3;
    t2=t2-t3;
  }
  /*pop out the stray field description if necessary*/
  if(stray_field) mcmagnet_pop(_particle);//this could be done with a stop bit in the magnet parameters
  
  /*inject another field description for the magnetic field betweeen pole shoes.*/
  mcmagnet_push(_particle,constant,&(ROT_A_CURRENT_COMP),&(POS_A_CURRENT_COMP),1,inside_magnet);
  
  /*see if we hit the magnet walls before we exit*/
  if (plane_intersect(&dt,x,y,z,vx,vy,vz,0,0,1,0,0,zdepth/2.0)==0){
    fprintf(stdout,"INFO: (%s) Missed the end plane of the magnet - this should not be\n",NAME_CURRENT_COMP);
    ABSORB;
  }
  if (fabs(dt-t4)>FLT_EPSILON) {
    if(verbose) fprintf(stdout,"INFO: (%s) hit magnet from inside gap.\n",NAME_CURRENT_COMP);
    mcmagnet_pop(_particle);
    ABSORB;
  }

  /*neutron is inside magnet gap - propagate to the foil plane or leave the neutron if it does not hit the plane
    while inside the gap*/
  if (foil_in){
    if (plane_intersect(&dt,x,y,z,vx,vy,vz,foil_n[0],foil_n[1],foil_n[2],0,0,0)==0 || dt<0 || dt>t4 ){
      if(verbose) fprintf(stderr,"INFO: (%s) Missed the foil flipper plane, (dt=%g).\n",NAME_CURRENT_COMP,dt);
      dt=dt;
      hit_foil=0;
    }
    else{
      PROP_DT(dt);
      t4=t4-dt;
      t2=t2-dt;
      /* at the foil - so flip spin/polarization*/
      if (foilthick == 0.0){
        foil_spin_flip(&sx,&sy,&sz,0.0,phi);
      }else{
        foil_spin_rot(&sx , &sy , &sz , &vx , &vy , &vz , phi , foilthick, By);
      }
      hit_foil=1;
    }
  }
  /*propagate to the other end of the magnet gap*/
  PROP_DT(t4);
  t2=t2-t4;
  /*pop the "inner" magnetic field*/
  mcmagnet_pop(_particle);
  
  if(stray_field)
    mcmagnet_push(_particle,0,&(ROT_A_CURRENT_COMP),&(POS_A_CURRENT_COMP),0,outstray_magnet);
  if (t2>1e-9){
    /*propagate to the end of the component*/
    PROP_DT(t2);
  }
  /*pop the 2nd stray field*/
  if(stray_field){
    mcmagnet_pop(_particle);
  }

%}

MCDISPLAY
%{
  double y=atan(phi)*zdepth/2.0;
  
  box(0.0,0.0,0.0,Bxwidth,Byheight,Bzdepth,0, 0, 1, 0);
  if (Bxwidth!=xwidth || Byheight!=yheight || Bzdepth!=zdepth){
    box(0,0,0,xwidth,yheight,zdepth,0, 0, 1, 0);
  }
  /*draw the foil plane*/
  multiline(5,-xwidth/2.0,-y,-zdepth/2.0, xwidth/2.0,-y,-zdepth/2.0, xwidth/2.0,y,zdepth/2.0, -xwidth/2.0,y,zdepth/2.0, -xwidth/2.0,-y,-zdepth/2.0);
%}
END