/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2011, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Instrument: Focalisation device
*
* %Identification
* Written by: <a href="mailto:desert@drecam.cea.fr">Sylvain Desert</a>
* Date: 2007.
* Origin: <a href="http://www-llb.cea.fr/">LLB</a>
* Release: McStas CVS_080624
* Version: $Revision: 1.4 $
* %INSTRUMENT_SITE: LLB
*
* Test instrument for neutron focalisation with a set of supermirrors. 
* No guide / velocity selector
* One parabolic SM converges the incoming beam to its focal point, then one 
* elliptic SM (with primary focal point at same location as parabolic SM focal point) 
* images the focal point onto the detector.
*
* %Description
* Test instrument for the MirrorElli and MirrorPara components
*
* Example: mcrun Focalisation.instr Qmin=0.0001 G=1 H=.0001 F=.00066 DET=8.2 lambda=14 divergence=.1 BeamWidth=.05 TetaMin=3
*
* %Parameters
* INPUT PARAMETERS:
* Qmin		: Minimum momentum transfer on the detector [AA^-1]
* G		: Ratio of the mirror focal length = increase in accepted divergence angle [1]
* H		: Height of mirrors [m]
* F		: Parabolic focal [m]
* DET 		: Distance between the two elliptic focal points [m]
* lambda   	: Mean wavelength of neutrons [AA]
* divergence	: source divergence (deg)
* BeamWidth	: Length of the source to focalize [m]
* TetaMin	: Minimum angle that must be reflected by the device assuming m=3 SM (TetaMin=3 for 5A and above neutron) [deg]
*
* %Link
* 
* %End
*******************************************************************************/
DEFINE INSTRUMENT Focalisation(Qmin=.0001,G=1,H=.0001,F=.00066,DET=8.2,lambda=14,divergence=.1,BeamWidth=.05,TetaMin=3)

DECLARE %{
  float T;       /* (Tan(2*Theta))^2 */
  float x1b;
  float y1b;
  float x2b;
  float y2b;
  float xf;
  float x1h;
  float y1h;
  float x2h;
  float y2h;
  float Ellia;	 /* Major axis of the ellipse */
  float Ellib;	 /* Minor axis of the ellipse */
  float BS;	 /* BeamStop size on the detector */
  float rayon1;	 /* Spatial filter radius */
  float A;
  float Delta;
%}

INITIALIZE %{
  T=(tan(TetaMin*3.1415926/180))*(tan(TetaMin*3.1415926/180));
  x1b=.5*(2*F+4*F/T+sqrt((2*F+4*F/T)*(2*F+4*F/T)-4*F*F));
  y1b=sqrt(4*F*x1b);
  y1h=y1b+BeamWidth;
  x1h=y1h*y1h/4/F;
  xf=(y1h*y1h-4*F*F)/4/F;
  Ellia=DET/2+G*F;
  Ellib=sqrt(Ellia*Ellia-(DET*DET/4));
  A=y1b/(x1b-F);
  Delta=(DET/Ellia/Ellia-2*A*A*F/Ellib/Ellib)*(DET/Ellia/Ellia-2*A*A*F/Ellib/Ellib)-4*(1/Ellia/Ellia+A*A/Ellib/Ellib)*(DET*DET/4/Ellia/Ellia+A*A*F*F/Ellib/Ellib-1);
  x2b=(2*A*A*F/Ellib/Ellib-DET/Ellia/Ellia-sqrt(Delta))/2/(1/Ellia/Ellia+A*A/Ellib/Ellib);
  A=y1h/(x1h-F);
  Delta=(DET/Ellia/Ellia-2*A*A*F/Ellib/Ellib)*(DET/Ellia/Ellia-2*A*A*F/Ellib/Ellib)-4*(1/Ellia/Ellia+A*A/Ellib/Ellib)*(DET*DET/4/Ellia/Ellia+A*A*F*F/Ellib/Ellib-1);
  x2h=(2*A*A*F/Ellib/Ellib-DET/Ellia/Ellia-sqrt(Delta))/2/(1/Ellia/Ellia+A*A/Ellib/Ellib);
  y2b=Ellib*sqrt((1-(x2b+Ellia)*(x2b+Ellia)/Ellia/Ellia));
  y2h=Ellib*sqrt((1-(x2h+Ellia)*(x2h+Ellia)/Ellia/Ellia));
  BS=Qmin*lambda/2/3.1415926*(DET+x2h);
  rayon1=BS*4*G*F/DET/T; 
%}

TRACE

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

COMPONENT source1 = Source_gen(dist = 0, xw = divergence , yh = .0001, Lambda0=lambda, 
			       dLambda=0.01, h = H/1000, w = BeamWidth) 
AT (0,0,0.1) RELATIVE PREVIOUS 

COMPONENT mirpara = MirrorPara(reflect="supermirror_m3.rfl", deltaxfaisceau = BeamWidth, ecartx=y1b,
				  yheight = H*2, focus = F) 
AT (-(y1b+BeamWidth/2), 0,0.1001) RELATIVE base

COMPONENT filter = Slit(xmin = 0, xmax = .1, ymin = -H, ymax = H)
AT ((rayon1+y1h)/(x1h-F)*(x1b-x1h)+y1h,0, x1h-x1b) RELATIVE PREVIOUS

COMPONENT spatialfilter = Slit(radius=rayon1) 
AT (0, 0, xf) RELATIVE mirpara

COMPONENT mirelli = MirrorElli(reflect="supermirror_m3.rfl",focus= F*G, interfocus = DET, yheight = H*2,
				  zmin=-x2b-DET/2, zmax=-x2h-DET/2)
AT (0, 0, DET/2-G*F+x1h) RELATIVE mirpara

COMPONENT monPSD = PSD_monitor(nx = 100, ny = 1, filename = "det.txt", xmin = -1.5*BS, xmax = 1.5*BS, 
			       yheight = .01)
AT (0, 0, DET/2) RELATIVE mirelli

FINALLY
%{
  printf("T=%f \n X1B=%f et Y1B=%f \n X1H=%f et Y1H=%f \n a=%f et b=%f \n X2H=%f et X2B=%f \n Y2H=%f et Y2B=%f\n\n",T,x1b,y1b,x1h,y1h,Ellia,Ellib,x2h,x2b,y2h,y2b);
  printf(" Delta=%f\n rayon1=%f\n",Delta,rayon1);
  printf(" Zmin=%f\n Zmax=%f\n",-x2b+x1h,-x2h+x1h);
  printf(" BS=%f\n xF=%f\n",BS,xf);
%}

END