S N Brown 2nd order achromat ;
;example 1 from K.L. Brown, SLAC-PUB-2257 (1979)
;in GICOSY H. Weick, 05.01.2013
; a focusing cell is defined with correctors
; after 4 such cells the system becomes a 2nd order achromat
; with a unity first order transfer matrix
;
 R P 9111.8 1 1 ; proton at 9.1 GeV adjusted to get correct Brho
 P X 0.03 0.003 ; initial phase space 30mm x 3 mrad in x
 P Y 0.03 0.003 ; initial phase space 30mm x 3 mrad in y
 D P 0.0 0.01 ;   deviation in momentum is 1% for plots
O C N M L N ;     output mode momentum and path length
;F T H ;
C M 1 ;           calculation mode fixed equations
;C M 2 ;           calculation mode differential algebra for fringe fields
;C M 3 ;           calculation mode differential algebra for all
;
C O 2 ;           calculation order 
O O 2 ;           output order 
;
S S ;             start system
;
; original values from K.L. Brown 1979
; MQX  =  8.02262 /10 ; kG -> T
; MQY  = -8.01363 /10 ; kG -> T
; MSX =  1.22642 /10 ; kG -> T
; MSY = -2.38768 /10 ; kG -> T
;
; refitted values in variables
 MQX =  8.014575556E-01 ;
 MQY = -8.021127910E-01 ;
 MSX =  1.226940521E-01 ;
 MSY = -2.390698123E-01 ;
;
; 
  ;F B MQX MQY ;   definition of a fit variables for 1st order fit
  F B MSX MSY ;    definition of a fit variables for 2nd order fit
  F D ;            disable fit
;
;-----------------------
B B EINS ;         block begin, named "EINS"
;
M Q 0.4 MQX 0.1 ;  magnetic quadrupole <length> <B-field> <radius>
M H 0.2 MSX 0.1 ;  magnetic hexapole <length> <B-field> <radius>
D L 0.3 ;          drift length
;
M S 166.782  2.0  0.1 ; magnetic sector field <radius> <angle> <half gap>
;
D L 0.3 ;
M H 0.2 MSY 0.1 ;
M Q 0.8 MQY 0.1 ;  
M H 0.2 MSY 0.1 ;
D L 0.3 ;
;
M S 166.782  2.0  0.1 ;
;
D L 0.3 ;
M H 0.2 MSX 0.1 ;  
M Q 0.4 MQX 0.1 ; 
;
B E ;              block end
;--------------------
;
B I EINS ;         apply block defined before
B I EINS ;
B I EINS ;
B I EINS ;
;
  F1 = ([X,X]-1.0)^2*9999 ;  construct variables for 1st order fit 
  F2 = ([Y,Y]-1.0)^2*9999 ;
  MIN1 = F1 + F2 ;           variable to minimize in fit
  ;P S MIN1,'|',MQX,' ',MQY ; print output to screen
  ;
  F4 = ([Y,BD]-0.0)^2*999 ;  construct variables for 2nd order fit 
  F5 = ([X,AD]-0.0)^2*999 ;
  MIN2 = F4 + F5 ;
  ;P S F4,' ',F5,' - ',MSX,' ',MSY ;
  F V MIN2 ;                 define fit variable to minimize
  F E 1E-20 1000 6 ;         fit <residuum> <interations> <fit method> 
;
D B 0.14 0.14 100 3 3 1 3 3 3 1 ;  draw beams
D S 0.14 10 100 3 3 1 3 1 1 ;      draw system with beams (top view)
D E 0.12 2.0 20 ([X,D]) ;	   draw expression
 D E 0.12 3000.0 20 ([X,AA]) ;
 D E 0.12 20.0 20 ([X,XA]) ;
 D E 0.12 20.0 20 ([X,AD]) ;
 D E 0.12 1.0 20 ([A,XD]) ;
 D E 0.12 200.0 20 ([Y,AY]) ;
 D E 0.12 50.0 20 ([Y,BD]) ;
 D E 0.12 1.0 20 ([B,YD]) ;
END ;                              end of calculation
;