INCLUDE 'COSY' ;
procedure run ;
   variable line 80 ;
   variable TRMAP_FT 1000 8 ;           {transport map from fp to target}
   variable TRMAP_TF 1000 8 ;           {transport map from target to fp}
   variable Q1MID_MAP 1000 8 ;           {temporary map to midpoint in Q1}
   variable Q2MID_MAP 1000 8 ;           {temporary map to midpoint in Q2}
   variable Q3MID_MAP 1000 8 ;           {temporary map to midpoint in Q3}
   variable TEMPMAP 1000 8 ;             {temporary map}
   variable RMAP 1000 4 ;                {reverse transport map 4 dim}
   variable TMPMAP 1000 8 ;
   variable TR1MAP 1000 8 ;
   variable MAPT 2000 5 ;
   variable OBJ 1 ;
   variable TMP1 1;
   variable TMPSPOS 1000;
   variable TDUMP_LEVEL 1;
   variable DIP_MOD 1 ;                 {dipole model number}

   variable t_ref 1 ;                   {kinetic energy of ref particle}
   variable p_cen 1 ;                   {central momentum setting in GeV}
   variable me 1 ;                       {mass of electron in GeV}
   variable mp 1 ;                       {mass of proton in GeV}
   variable mid 1 ;
   variable exit 1 ;

   variable RES 1 4 ; 
   variable RORDER 1 ;                   {order of recon matrix}
   variable FORDER 1 ;                   {order of forward matrix}
   variable FRMODE 1 ;                   {fringe field mode}

{resolutions for reverse map}
   variable R1 1 ;                       {x resolution}
   variable R2 1 ;                       {dx/dz resolution}
   variable R3 1 ;                       {y resolution}
   variable R4 1 ;                       {dy/dz resolution}
   variable R5 1 ;                       {delta resolution}
   variable R6 1 ;                       {detector position resolution}
   variable R7 1 ;                       {detector angle resolution}

{multipole variables}
   variable Q1 1 ;            { Quadrupole strength Q1 }
   variable H1 1 ;
   variable O1 1 ;
   variable D1 1 ;
   variable DD1 1 ;
   variable Q2 1 ;            { Quadrupole strength Q2 }
   variable H2 1 ;
   variable O2 1 ;
   variable D2 1 ;
   variable DD2 1 ;
   variable Q3 1 ;            { Quadrupole strength Q3 }
   variable H3 1 ;
   variable O3 1 ;
   variable D3 1 ;
   variable DD3 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_N 1 5 ;       { field index area }
   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 ;


{procedure declarations}

{-----------------------------------------------------------------------------}

   procedure PMT IU ;
   variable I 1 ;
      LOOP I 1 5 ;
       MAPT(I) := TMPMAP(I) ;
      ENDLOOP ;
      DAPRV Mapt 5 5 TWOND IU ;
   ENDprocedure ;

{-----------------------------------------------------------------------------}

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

{-----------------------------------------------------------------------------}

   procedure SET_ORDERS ;
     RORDER:=5 ;
     FORDER:=5 ;
     FRMODE:=3 ;
   endprocedure ;

{-----------------------------------------------------------------------------}

   procedure DEF_QUADS IU ;
{      Q1 := -1.16452 ; }
{      Q1 := -1.16948 ; }
      Q1 := -0.243838 ;
      H1 := 0.0 ; 
      O1 := 0.0 ; 
      D1 := 0.0 ; 
      DD1 := 0.0 ;

{      Q2 := 0.923703 ; }
{      Q2 := 0.92556 ; }
      Q2 := 0.1934 ;
      H2 := 0.0 ; 
      O2 := 0.0 ; 
      D2 := 0.0 ; 
      DD2 := 0.0 ;

{      Q3 := 0.854484 ; }
{      Q3 := 0.89647 ; }
      Q3 := 0.17892 ;
      H3 := 0.0 ; 
      O3 := 0.0 ;
      D3 := 0.0 ; 
      DD3 := 0.00245 ;

      IF IU=3 ;
         FC 2 1 1 0.1122 8.5000 -1.4982 3.5882 -2.1209 1.723 ;
         FC 2 2 1 0.1122 8.5000 -1.4982 3.5882 -2.1209 1.723 ;
         FC 3 1 1 0.1122 8.5000 -1.4982 3.5882 -2.1209 1.723 ;
         FC 3 2 1 0.1122 8.5000 -1.4982 3.5882 -2.1209 1.723 ;
         FC 4 1 1 0.1122 8.5000 -1.4982 3.5882 -2.1209 1.723 ;
         FC 4 2 1 0.1122 8.5000 -1.4982 3.5882 -2.1209 1.723 ;
         FC 5 1 1 0.1122 8.5000 -1.4982 3.5882 -2.1209 1.723 ;
         FC 5 2 1 0.1122 8.5000 -1.4982 3.5882 -2.1209 1.723 ;
         FC 6 1 1 0.1122 8.5000 -1.4982 3.5882 -2.1209 1.723 ;
         FC 6 2 1 0.1122 8.5000 -1.4982 3.5882 -2.1209 1.723 ; 
      ELSEIF 1=1 ;
         FC 2 1 1 0.1122 6.2671 -1.4982 3.5882 -2.1209 1.723 ;
         FC 2 2 1 0.1122 6.2671 -1.4982 3.5882 -2.1209 1.723 ;
         FC 3 1 1 0.1122 6.2671 -1.4982 3.5882 -2.1209 1.723 ;
         FC 3 2 1 0.1122 6.2671 -1.4982 3.5882 -2.1209 1.723 ;
         FC 4 1 1 0.1122 6.2671 -1.4982 3.5882 -2.1209 1.723 ;
         FC 4 2 1 0.1122 6.2671 -1.4982 3.5882 -2.1209 1.723 ;
         FC 5 1 1 0.1122 6.2671 -1.4982 3.5882 -2.1209 1.723 ;
         FC 5 2 1 0.1122 6.2671 -1.4982 3.5882 -2.1209 1.723 ;
         FC 6 1 1 0.1122 6.2671 -1.4982 3.5882 -2.1209 1.723 ;
         FC 6 2 1 0.1122 6.2671 -1.4982 3.5882 -2.1209 1.723 ; 
      ENDIF ;
   ENDprocedure ;

{-----------------------------------------------------------------------------}

   procedure DEF_DIPOLE ;
   variable I 1 ;
   IF DIP_MOD<0.5 ;
      D_RADIUS := 8.40 ; D_BEND := 45.0 ; D_GAP := 0.125 ;
      D_ALPHA := -30.0 ; D_BETA := -30.0 ;
      D_N(1) := -1.25 ;
      D_N(2) := -D_N(1)*D_N(1) ; 
      D_N(3) := D_N(1)*D_N(1)*D_N(1) ;
      D_N(4) := -D_N(1)*D_N(1)*D_N(1)*D_N(1) ; 
      D_N(5) := D_N(1)*D_N(1)*D_N(1)*D_N(1)*D_N(1) ;
      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 ;
   ELSEIF 1>0 ;
      D_RADIUS:=8.458 ;
      D_BEND:=45.0 ;
      D_GAP:=0.1244 ;
      D_H1:=0.6369 ; D_H2:=-0.3571 ;
      D_ALPHA := -30.0 ; D_BETA := -30.0 ;
      D_N(1) := -1.26 ; 
      D_N(2) := -D_N(1)*D_N(1) ; 
      D_N(3) := D_N(1)*D_N(1)*D_N(1) ;
      D_N(4) := -D_N(1)*D_N(1)*D_N(1)*D_N(1) ; 
      D_N(5) := D_N(1)*D_N(1)*D_N(1)*D_N(1)*D_N(1) ;
      D_S1(1) := -0.5513012 ; D_S1(2) := 0.4656544 ; 
      D_S1(3) := 0.03799517 ; D_S1(4) := 0.03613541 ; 
      D_S1(5) := -0.2209376;
      D_S2(1) := -0.6049777 ; D_S2(2) := 0.2850934 ;
      D_S2(3) := 0.4569041 ; D_S2(4) := 0.4247217 ;
      D_S2(5) := -0.9568408 ;
   ENDIF ;

   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 ;

   endprocedure ;

{********************************************************************}


   OV 5 3 0 ;                          {Order Value <order><phase dim><#par>}
{   RP 4000.0 0.0005 -1;}
   RPE 837.56 ;
   PTY 0 ;
   SB .0 .030 .0 .05 .020 .0 .0 .05 .0 .0 .0 ;

   SET_ORDERS ;
   TDUMP_LEVEL:=3;
   DIP_MOD:=0;
   OPENF 20 'hrs_forward_map.meek' 'UNKNOWN' ; 

   DEF_QUADS 1 ; DEF_DIPOLE ;
   UM ;

{-------midpoint maps in Q1 and Q2--------------}

  IF TDUMP_LEVEL=3 ;
   DL 1.6000 ;
   FR 1 ;
   mid:=0.9413*2.0/3.0 ;
   MQ mid Q1 0.1500 ;
   SM Q1MID_MAP ;
   exit:=0.9413*0.3333333;
   MQ exit Q1 0.1500 ;
   
   DL 1.25 ;

   mid:=1.8266*0.666667 ;
   M5 mid Q2 H2 O2 D2 DD2 0.300 ; {Q2 2/3}
   SM Q2MID_MAP ;                              {transport map at 2/3 of Q2}
  ENDIF ;

{-----------------------------------------------}
 
      FR FRMODE ;
      UM ; CR ; ER 1 3 2 3 1 3 1 1 ; BP ;

      DL 1.5903 ;                         {Drift to Q1}

     DEF_QUADS 0 ;
     IF TDUMP_LEVEL>1 ;
      PT 20 ;
     ENDIF ;

     IF TDUMP_LEVEL=3 ;                  {PT midpoint of Q1}
      SM TMPMAP ;
      UM ;
      AM Q1MID_MAP ;
      PT 20 ;
      UM ;
      AM TMPMAP ;
     ENDIF ;
     FR -1 ; MQ 0.9413 Q1 0.1500 ;      {Q1 entrance fringe}
     FR  0 ; MQ 0.9413 Q1 0.1500 ;      {Q1 central}
     FR -2 ; MQ 0.9413 Q1 0.1500 ;      {Q1 exit fringe}

     IF TDUMP_LEVEL>1 ;
      PT 20 ;
     ENDIF ;

      DL 1.172 ;                          {Drift to Q2}

     DEF_QUADS 0 ;
     IF TDUMP_LEVEL>1 ;
      PT 20 ;
     ENDIF ;

     IF TDUMP_LEVEL=3 ;                     {PT midpoint of Q2}
      SM TMPMAP ;
      UM ;
      AM Q2MID_MAP;
      PT 20 ;
      UM ;
      AM TMPMAP ;
     ENDIF ;


     DEF_QUADS 0 ;
     FR -1 ; M5 1.8266 Q2 H2 O2 D2 DD2 0.300 ; {Q2 entrance fringe}
     FR  0 ; M5 1.8266 Q2 H2 O2 D2 DD2 0.300 ;   {Q2 central}
     FR -2 ; M5 1.8266 Q2 H2 O2 D2 DD2 0.300 ; {Q2 exit fringe}
{      M5 1.8266 Q2 H2 O2 D2 DD2 0.300 ;  }   {Q2}

     IF TDUMP_LEVEL>1 ;
      PT 20 ;
     ENDIF ;

      DL 4.4308 ;                          {Drift to B1}

     IF TDUMP_LEVEL>1 ;
      PT 20 ;
     ENDIF ;

     DEF_DIPOLE ;
     FR FRMODE ;
     MC D_RADIUS D_BEND D_GAP D_N D_S1 D_S2 5 ;          {B1}

{      MC D_RADIUS D_BEND D_GAP D_N D_S1 D_S2 5 ;} {B1}

     IF TDUMP_LEVEL>1 ;
      PT 20 ;
     ENDIF ;

      DL 1.5925 ;                         {Drift to Q3}

     IF TDUMP_LEVEL>1 ;
      PT 20 ;
     ENDIF ;

     DEF_QUADS 0 ;
     FR -1 ; M5 1.8268 Q3 H3 O3 D3 DD3 0.300 ; {Q3 entrance fringe}
     FR  0 ; M5 1.8268 Q3 H3 O3 D3 DD3 0.300 ;   {Q3 central}
     FR -2 ; M5 1.8268 Q3 H3 O3 D3 DD3 0.300 ; {Q3 exit fringe}

{      M5 1.8268 Q3 H3 O3 D3 DD3 0.300 ;} {Q3}

     IF TDUMP_LEVEL>1 ;
      PT 20 ;
     ENDIF ;

      DL 3.4523 ;                          {Drift to detector plane}

      PS .50 ;                             {Mark focal plane}
      DL 5.000 ;                           {extend rays for picture}
      PS .50 ;                             {Mark plane}
      EP ;
      DL -5.000 ;                          {Drift back to focal plane}
      EP ;

     IF TDUMP_LEVEL>1 ;
      PT 20 ;
     ENDIF ;
 
   CLOSEF 20 ;

   OPENF 20 'hrs_reconstruction_map.meek' 'UNKNOWN' ; 

   TMPSPOS := SPOS ; 
   CATR MAP TR1MAP SPOS ; 
   UM ; 
   SPOS := TMPSPOS ; 
   AM TR1MAP ;

   RR MAP 0 0.01 0.01 0.01 0.01 0.0 0.0 1000 RORDER RMAP RES ;

   pmr 20 ;
   closef 20 ;

   PG -10 -10 ;
   PP -61 0 90 ;
   PP -62 0 0 ;
endprocedure ;

run ;

end ;