{------------------------------------------------------------------------------
 HRS-MONTE.FOX

 Generate Monte-Carlo forward and reconstruction matrix elements.

 Hall A HRS spectrometer.

 All output matrix elements are in TRANSPORT coordinates.
 Note that COSY-7 uses different units than COSY-6 when printing the maps
 in transport notation (procedure PT.) The reconstruction map is still
 calculated using the old COSY-6 units. The output format corresponds to what
 is expected by my program monte_hms.

 X-dependent (raster) reconstruction map generated.

 Version: 0.0-dev
 Author: D. Potterveld, Argonne National Lab, 26-Jun-2001

 Modification History:

 ------------------------------------------------------------------------------
}
INCLUDE 'COSY' ;
PROCEDURE RUN ;

{=========================== Variable Declarations ===========================}

   VARIABLE Z1C 1 ; VARIABLE L1 1 ; VARIABLE L1E 1 ;
   VARIABLE Z2C 1 ; VARIABLE L2 1 ; VARIABLE L2E 1 ;
   VARIABLE Z3C 1 ; VARIABLE L3 1 ; VARIABLE L3E 1 ;
   VARIABLE ZDI 1 ;
   VARIABLE DR1 1 ;
   VARIABLE DR2 1 ;
   VARIABLE DR3 1 ;
   VARIABLE DRD 1 ;
   VARIABLE Q1 1 ;        { quadrupole strength }
   VARIABLE Q2 1 ;        { quadrupole strength }
   VARIABLE Q3 1 ; VARIABLE DD3 1 ; {Q3 parameters }

   VARIABLE Q1MID 1 ; VARIABLE Q2MID 1 ; VARIABLE Q3MID 1 ;

{dipole variables}

   variable D_RADIUS 1 ;    { dipole bend radius }
   variable D_BEND 1 ;      { dipole bend angle }
   variable D_GAP 1 ;       { dipole half gap }
   variable D_ALPHA 1 ;     { entrance shim angle }
   variable D_BETA 1 ;      { exit shim angle }
   variable D_INDEX 1 ;     { Field index }
   variable D_N 1 5 ;       { field radial dependence coefficients }
   variable DELTA_D_N 1 ;   { error on the field index }
   variable D_S1 1 5 ;      { entrance polynomial coefficients }
   variable D_S2 1 5 ;      { exit polynomial coefficients }
   variable D_H1 1 ;        { dipole entrance curvature }
   variable D_H2 1 ;        { dipole exit curvature }
   variable D_S01 1 5 ;
   variable D_S2 1 5 ;
   variable D_S02 1 5 ;

   VARIABLE OBJ 1 ;       { function used in evaluating first order focus }
   VARIABLE FILE 80 ;     { filename }
   VARIABLE LINE 80 ;     { dummy for user to answer into }
   VARIABLE TILT 1 ;      { Focal plane tilt angle }
   VARIABLE I 1 ;
   VARIABLE J 1 ;
   VARIABLE TMP 1 ;
   VARIABLE SNAME 3 ;     {Spectrometer name}
   VARIABLE FRMODE 1 ;    { Fringe field mode }
   VARIABLE MORDER 1 ;    { Maximum order of maps }
   VARIABLE FORDER 1 ;    { Order of forward maps }
   VARIABLE RORDER 1 ;    { Order of reconstruction map }
   VARIABLE RES 1 4 ;     { Resolution array }

   VARIABLE TMPSPOS 1000 ;  { Saved value of SPOS }
   VARIABLE POSITION 1 ;    { Position of each output plane }
   VARIABLE GMAP 2000 8 ;   { Global transfer map }
   VARIABLE RMAP 2000 5 ;   { Reconstruction map }

{========================== Procedure declarations ===========================}

 { Print rmap to IU }

   PROCEDURE PMR IU ; 
      VARIABLE I 1 ; VARIABLE M 2000 5 ; 
      LOOP I 1 5 ; M(I) := RMAP(I) + 0*DD(1) ; ENDLOOP ;
      DAPRV M 4 5 TWOND IU ; ENDPROCEDURE ;

 { Reconstruction with X as independent variable }

   PROCEDURE RR1 M X0 X A Y B D PR AR N OR MR RES ;{TRAJECTORY RECONSTRUCTION}
      PROCEDURE REC ; VARIABLE NN 1 ; VARIABLE S MAX(NM1,N+1) ;
      VARIABLE I 1 ; VARIABLE M4 NM1 5 ; VARIABLE CM 1 6 ;
      VARIABLE RI N 6 ; VARIABLE RO N 4 ; VARIABLE RB N 4 ; VARIABLE RE N 4 ;
      VARIABLE IRAN 1 ; VARIABLE J NM1 ; VARIABLE NOM 1 ; VARIABLE IER 1 ;
      VARIABLE MR4 NM1 4 ;
      FUNCTION RAN IRAN ; RAN := 98765*IRAN+.12345 ; RAN := 2*(RAN-NINT(RAN)) ;
         IRAN := RAN ; ENDFUNCTION ; IRAN := 0 ;
      IF OR>0 ; LOOP I 1 4 ; CM(I+1) := DD(I) ;
         ENDLOOP ; NOM := NOC ; CO OR ; CM(1) := DD(5) ; CM(6) := 0*DD(1) ;
         IF ND=3 ; CM(6) := DD(4) ; CM(5) := 0*DD(1) ; ENDIF ;
         POLVAL 1 M 4 CM 6 M4 4 ; 
         M4(5) := DD(5) ; MI M4 MR 5 IER ; CO NOM ; ENDIF ;
      S := 0 ; LOOP I 1 N-1 ; S := S&0 ; ENDLOOP ;
      LOOP I 1 N ; VELSET S I X0+X*RAN(IRAN) ; ENDLOOP ; RI(1) := S ;
      LOOP I 1 N ; VELSET S I A*RAN(IRAN) ; ENDLOOP ; RI(2) := S ;
      LOOP I 1 N ; VELSET S I Y*RAN(IRAN) ; ENDLOOP ; RI(3) := S ;
      LOOP I 1 N ; VELSET S I B*RAN(IRAN) ; ENDLOOP ; RI(4) := S ;
      LOOP I 1 N ; VELSET S I D*RAN(IRAN) ; ENDLOOP ; RI(5) := S ;
      RI(6) := S ;  POLVAL 1 M 4 RI 6 RO 4 ;
      LOOP I 1 N ; VELSET S I PR*RAN(IRAN) ; ENDLOOP ; RE(1) := S ;
      LOOP I 1 N ; VELSET S I AR*RAN(IRAN) ; ENDLOOP ; RE(2) := S ;
      LOOP I 1 N ; VELSET S I PR*RAN(IRAN) ; ENDLOOP ; RE(3) := S ;
      LOOP I 1 N ; VELSET S I AR*RAN(IRAN) ; ENDLOOP ; RE(4) := S ;
      RO(1) := RO(1) + RE(1) ; RO(2) := RO(2) + RE(2) ;
      RO(3) := RO(3) + RE(3) ; RO(4) := RO(4) + RE(4) ;
      LOOP I 1 4 ; DAPLU MR(I) 5 X0 S ; MR4(I) := S ; ENDLOOP ;
      POLVAL 1 MR4 4 RO 4 RB 4 ;  RE(1) := RI(2)-RB(1) ; RE(2) := RI(3)-RB(2) ;
      RE(3) := RI(4)-RB(3)  ; RE(4) := RI(6)-RB(4) ;
      LOOP J 1 4 ; RES(J) := 0 ; NN := 0 ; LOOP I 1 N ; VELGET RE(J) I S ;
         IF S#0 ; RES(J) := RES(J) + ABS(S) ; NN := NN + 1 ; ENDIF ; ENDLOOP ;
         RES(J) := NN/(RES(J)+1E-20)/4 ; IF IER#0 ; RES(J) := 0 ; ENDIF ;
         ENDLOOP ; ENDPROCEDURE ;
      NM1 := MAX(NM1,N) ; SCRLEN 2*NM1 ; REC ;
      NM1 := NMON(NO+1,NV) ; SCRLEN LSCR ; ENDPROCEDURE ;

 { Ask user for input value }

   PROCEDURE ASK VAR PROMPT LL UL ;
      VARIABLE TMP 1 ;
      VAR := 2*ABS(UL) + 1 ;
      WHILE (VAR<LL)+(VAR>UL) ; WRITE 6 '$'&PROMPT ; READ 5 VAR ;
      ENDWHILE ; ENDPROCEDURE ;

 { Output forward transformation data at an aperture and update internal maps }

   PROCEDURE APERTURE LABEL REGION ;
      WRITE 20 'NAME: '&LABEL ;
      WRITE 20 'REGION: '&REGION ;
      WRITE 20 'OFFSET: '&SF(TMPSPOS,'(F25.16)')&' (in meters)' ;
      WRITE 20 'LENGTH: '&SF(SPOS,'(F25.16)')&' (canonical length in meters)' ;
      PT 20 ;                           {Output data}
      TMPSPOS := TMPSPOS + SPOS ;
      ANM MAP GMAP GMAP ; UM ;          {Update global map + reset}
      WRITE 6 ' Aperture at S = '&S(TMPSPOS) ;
   ENDPROCEDURE ;

 { Compute and output maps to inner aperture and exit of HRS quad magnet. }

   PROCEDURE HRSQUAD NAME LEFF Q HEX OCT DEC DDEC A FRAC ;
      VARIABLE DRIFT 1 ; VARIABLE LMID 1 ; VARIABLE QFLAG 1 ;

{ Compute flag for element type }

      QFLAG := ABS(HEX)+ABS(OCT)+ABS(DEC)+ABS(DDEC);
      IF QFLAG=0 ;  WRITE 6 NAME&': Using normal quad' ;
      ELSEIF TRUE ; WRITE 6 NAME&': Using M5 quad' ;
      ENDIF ;

{ Compute the location of the inner aperture }

      LMID  := (LEFF*FRAC) ;

      IF FRMODE#0 ; FR -1 ;                       {Entrance fringe field}
        IF QFLAG=0 ;  MQ LEFF Q A ;
        ELSEIF TRUE ; M5 LEFF Q HEX OCT DEC DDEC A ; ENDIF ;
      ENDIF;

      FR  0 ;                                     {First part of quad}
      IF QFLAG=0 ;  MQ LMID Q A ;
      ELSEIF TRUE ; M5 LMID Q HEX OCT DEC DDEC A ; ENDIF ;

      APERTURE NAME&'_MID' 'FRONT' ;

      IF QFLAG=0 ;  MQ (LEFF-LMID) Q A ;            {Second part of quad}
      ELSEIF TRUE ; M5 LEFF-LMID Q HEX OCT DEC DDEC A ; ENDIF ;

      IF FRMODE#0 ; FR -2 ;                       {Exit fringe field}
        IF QFLAG=0 ;  MQ LEFF Q A ;
        ELSEIF TRUE ; M5 LEFF Q HEX OCT DEC DDEC A ; ENDIF ;
      ENDIF;

      APERTURE NAME&'_EXIT' 'BACK' ;

      FR FRMODE ; ENDPROCEDURE ;                  {Restore fringe field mode}


{================================= Main Code =================================}

{ Computation order specifications }

   MORDER := 5 ;                        {Max order for calculations}
   ASK FORDER 'Order for maps (1-'&SF(morder,'(i1)')&'): ' 1 morder ;

   OV FORDER 3 0 ;                      {Order Value <order><phase dim><#par>}
   CO FORDER ; RORDER := FORDER ;

{ Beam specifications }

   RPE 837.56 ;                          {REF PARTICLE <Momentum>}

   SB .01 .030  0   .05 .020  0   0  .05  0   0  0 ;
{  --> PX   PA  R12  PY   PB  R34 PT  PD  R56 PG PZ      }


{ Spectrometer name }

   SNAME := 'HRS' ;

{ Fringe-field specifications.}

   ASK frmode 'Fringe-field mode (0-3): ' 0 3 ; 
   FR FRMODE ;

{ From Dave Meekins HRS file }

   FC 1 1 1 0.04725 2.2395 -0.9768 0.7288 -0.1299 0.0222 ;
   FC 1 2 1 0.04725 2.2395 -0.9768 0.7288 -0.1299 0.0222 ;

   LOOP I 2 6 ; LOOP J 1 2 ;
      FC I J 1 0.1122 6.2671 -1.4982 3.5882 -2.1209 1.723 ;
   ENDLOOP ; ENDLOOP ;

{ Magnet geometry}

   L1 := 0.9413 ;            {Magnetic length of Q1}
   L2 := 1.8266 ;            {Magnetic length of Q2}
   L3 := 1.8268 ;            {Magnetic length of Q3}

   DR1 := 1.5903 ;                  {Drift from target to Q1}
   DR2 := 1.172  ;                  {Drift from Q1 to Q2}
   DRD := 4.4308 ;                  {Drift from Q2 to dipole}
   DR3 := 1.5925 ;                  {DRIFT FROM DIPOLE TO Q3}

{ QUAD STRENGTHS }

   Q1 := -0.243838 ;
   Q2 :=  0.1934 ; 
   Q3 :=  0.17892 ;
   DD3 := 0.00245 ;

{ Dipole parameters }

      D_RADIUS := 8.40 ; D_BEND := 45.0 ; D_GAP := 0.125 ;
      D_ALPHA := -30.0 ; D_BETA := -30.0 ; D_INDEX := -1.25 ;

      D_N(1) := D_INDEX ;
      loop i 2 5 ; d_n(i) := -d_n(i-1)*D_INDEX ; endloop ;

      D_S1(1) := TAN(-30.*DEGRAD) ;   {Entrance pole shape}
      D_S1(2) :=  0.4014526 ;
      D_S1(3) :=  0.21710 ;
      D_S1(4) := -0.0825 ;
      D_S1(5) := -0.3110 ;
      D_S2(1) := TAN(-30.*DEGRAD) ;   {Exit pole shape}
      D_S2(2) := -0.4433318 ;
      D_S2(3) :=  0.18020 ;
      D_S2(4) := -0.2443 ;
      D_S2(5) :=  1.1649 ;


{ Begin composing forward maps }

   file := SNAME&'_fr'&sf(frmode,'(i1)')&'_for_maps.dat' ;
   openf 20 file 'NEW' ;

   FR FRMODE ;
   UM ; SM GMAP ; TMPSPOS := SPOS ;  {Initialize global map}

{ FRONT QUADS }

      DL DR1 ;                   {Drift to Q1 physical aperture}
      APERTURE 'Q1_ENTRANCE' 'FULL' ;
      HRSQUAD  'Q1' L1 Q1 0 0 0 0 0.15 2/3 ;   {Q1}

      DL DR2 ;                        {Drift to Q2}
      APERTURE 'Q2_ENTRANCE' 'FULL' ;
      HRSQUAD  'Q2' L2 Q2 0 0 0 0 0.30 2/3 ;   {Q2}

{ DIPOLE }

      DL DRD ;                        {Drift to dipole}
      APERTURE 'DIPOLE_ENTRANCE' 'FULL' ;

{ Consistency check: we better be at the right Z! }

{
      IF ABS(TMPSPOS-ZDI)>1E-9 ;
        WRITE 6 'Warning: Dipole location is bad:' ;
        WRITE 6 'POSITION = '&S(TMPSPOS)&', but should be '&S(ZDI) ;
      ENDIF ;
}
      MC D_RADIUS D_BEND D_GAP D_N D_S1 D_S2 5 ;          {B1}

      APERTURE 'DIPOLE_EXIT' 'FULL' ;

{ REAR QUAD }

      DL DR3 ;                   {Drift to Q1 physical aperture}
      APERTURE 'Q3_ENTRANCE' 'FULL' ;
      HRSQUAD  'Q3' L3 Q3 0 0 0 DD3 0.30 2/3 ;   {Q3}

{ Reconstruction plane }

      DL 3.4523 ;                        {Drift to detector plane}
      APERTURE 'RECON PLANE' 'FULL' ;

   closef 20 ;                          {Done outputing forward MAPS}
{ Convert GMAP from canonical to TRANSPORT form}

   CATR GMAP MAP TMPSPOS ; spos := tmpspos ;

{ Compute focal objective function }

   OBJ := SQRT(ME(1,2)*ME(1,2)+ME(3,4)*ME(3,4)) ;
   WRITE 6 ' Q1,Q2, 1st order OBJ = ' Q1&Q2&OBJ ;
   WRITE 6 ' Total length = '&S(CONS(SPOS))&' meters' ;

{ Compute focal plane tilt angle }

   TILT := (180*atan(ME(1,1)*ME(1,26)/ME(1,6))/pi)+90 ;
   write 6 ' Focal plane tilt ='&SF(tilt,'(f8.3)')&' degrees' ;

{ Compute RECON map }

   RR1 MAP 0 0 .04 0 .04 .15 0 0 1000 RORDER RMAP RES ;
   write 6 ' RECON resolutions =  ' RES(1)&RES(2)&RES(3)&RES(4) ;
   file := SNAME&'_fr'&sf(frmode,'(i1)')&'_rec.dat' ;
   openf 20 file 'NEW' ; pmr 20 ; closef 20 ;

ENDPROCEDURE ;
RUN ;
END ;