/ ------------------------------------------------------------------
/ ---------        Model for CLEO Solenoidal Field         ---------
/ ---------                                                ---------
/ ---------               J.A.Crittenden                   ---------
/ ---------                 12 May 2004                    ---------
/ ------------------------------------------------------------------
/
/ 12.05.2004  Version 58   First 2D model of CLEO solenoidal field. Simple pole tip. Half coil. No symmetry around X axis.
/ 12.05.2005  Version 59   Use POLYGON to make air/steel in order to choose element size. Simple pole. Beautiful model.
/ 12.05.2005  Version 60   Put in CLEO II and CLEO III pole options. Choose CLEO II. 3256A for comparison with IEEE paper.
/                          This version gives mesh errors for the CLEO III pole tip geometry.
/ 13.05.2005  Version 61   Modified pole tip air volumes for CLEO III. Added cleo version and current to file names.
/ 13.05.2005  Version 62   Corrected mistake in bias of background side at negative Y. 
/                          Decreased mesh size around pole tips and in beam region.
/ 
/ 24.01.2011  Version 63   Update to Version 13.0
/
/ ============================================================ Parameters ===========================
/                                                        Units are cm, amps
/ ===================================================================================================
UNITS LENGTH=CM FLUX=GAUSS FIELD=OERSTED POTENTIAL=GCM CONDUCTIVITY=SCM DENSITY=ACM2 FORCE=DYNE ENERGY=ERG POWER=ERGS MASS=GRAMME

/ ---- Model version number ----
$ string version '63'

/ --- Solenoid operating current ---
/$ string current '3256'
$ string current '2200'
$ const #solcurrent -&current&

/ CLEO pole tip version
$const #cleo3 1
$if #cleo3 eq 1
$ string cleover 'III'
$else 
$ string cleover 'II'
$end if

/ --- Choose field calculation method ---
$const #integrate 1

/ --- Enable/disable table file writing ---
$const #table 1

/ --- Choose mesh scale size ---  A value of 1 chooses 10-cm element size
$ const #meshscale 1.0

/ --- Choose interactive/batch mode ---
$ const #interact 1


$if #integrate eq 1
/ To choose field calculation by integration over currents and magnetizations, use this line
 SET FIELD=MAGNETIC +INTEGRAL FLDTYPE=TOTAL SYMMETRY=AXISYMMETRY SOLUTION=AXI  ELEMENT=LINEAR
$else
 SET FIELD=MAGNETIC  SYMMETRY=AXISYMMETRY SOLUTION=AXI  REFXAXIS=INVE ELEMENT=LINEAR
$end if

/
/ The CLEO solenoid has 1300 turns in two radial layers. They are arranged longitudinally such that
/ the wires in the central region are spaced by 5.1 mm and those in the two end regions spaced by 4.89 mm.
/ The only information I have found on the lengths of these regions is in an unpublished, undated note
/ which states "The winding density is about 4% great on the outer quarters ..."
/ One can calculate the lengths via (x1/0.51 + x2/0.489) = 650/2 turns together with a ratio x1/x2 which
/ gives the desired field uniformity of 10^-3 over the DR volume. Using a value of x1/x2=1.08 taken 
/ from Oxford publication "Design Study of a Large superconducing Magnet CLEO 2 for Cornell University" (5 Oct 84)
/ we obtain x1=845 mm, x2=780 mm. This publication also provided the B-H curve for the CLEO yoke.
/ But this calculation gives a total length of only 3.2 m, and all sources give the active length as 3.5 m.
/ Also, the field measurements in the IEEE paper cannot be reproduced with a coil length of 3.2 m, even when
/ the distance of the pole tip is varied from 1.61 to 1.81 m (studies done with version 15).
/ So with version 16, we set the coil length to 3.5 m, the density change at 1/4 length, and allow
/ the two current densities to vary from those corresponding to 5.1 and 4.89 mm. They become something
/ more like 5.3 mm spacing.
/
/ Further parameters of the coil were obtained from IEEE Trans Act Nucl Sci Vol 37, No 3, p. 1172 (1990)
/ So were the barrel yoke plate thicknesses.
/
/ ---- Solenoid parameters ----
/ Coil thickness
$ const #soldr 2.0
/ Coil inner radius
$ const #solradin 154.
$ const #solradout #solradin+#soldr
/ Change the shape of this field with this parameter, which is the half-length of the region with lower current density.
/ A smaller value raises the field further from the center of the solenoid.
/ The calculation keeps the total number of turns at 1300 in two radial layers while changing the total length of the coil.
$ const #solmidyhi 86.
$ const #solmidylo -#solmidyhi
$ const #enddens 0.489
$ const #middens 0.510
$ const #midturns 1300./(1+#middens/#enddens)
/ This factor of two is for the two end coil regions
$ const #endturns (1300-#midturns)/2
$ const #solendyhi 175
$ const #solendylo -#solendyhi
$ const #solmiddy #solmidyhi-#solmidylo
$ const #midcurrd #solcurrent*#midturns/(#soldr*#solmiddy)
$ const #solenddy #solendyhi-#solmidyhi
$ const #endcurrd #solcurrent*#endturns/(#soldr*#solenddy)
/
/============================================================================================================================
/ Create three solenoids of differing current density 
DRAW +DEFAULTS SHAPE=POLYGON TOLERANCE=5.0E-05 MATERIAL=1 PERM=1 DENS=#midcurrd CONDUCTIVITY=0 PHASE=0 VELOCITY=0 N=0 SYMMETRY=0 XCEN=0 YCEN=0 ANGLE=0 NX=1 DX=0 NY=1 DY=0 ROTATIONS=1 TROTATION=0 MIRROR=NO
DRAW X12=#solradin Y14=-#solmidyhi X34=#solradout Y23=#solmidyhi shape=Q N1=10*#meshscale N2=2*#meshscale N3=10*#meshscale N4=2*#meshscale B1=0.5 B2=0.5 B3=0.5 B4=0.5 

DRAW +DEFAULTS SHAPE=POLYGON TOLERANCE=5.0E-05 MATERIAL=1 PERM=1 DENS=#endcurrd CONDUCTIVITY=0 PHASE=0 VELOCITY=0 N=0 SYMMETRY=0 XCEN=0 YCEN=0 ANGLE=0 NX=1 DX=0 NY=1 DY=0 ROTATIONS=1 TROTATION=0 MIRROR=NO
/DRAW X12=#solradin Y14=#solendylo X34=#solradout Y23=#solmidylo shape=Q N1=10*#meshscale N2=2*#meshscale N3=10*#meshscale N4=2*#meshscale B1=0.5 B2=0.5 B3=0.5 B4=0.5 
DRAW X12=#solradin Y14=#solmidyhi X34=#solradout Y23=#solendyhi shape=Q N1=10*#meshscale N2=2*#meshscale N3=10*#meshscale N4=2*#meshscale B1=0.5 B2=0.5 B3=0.5 B4=0.5 
DRAW X12=#solradin Y14=-#solendyhi X34=#solradout Y23=-#solmidyhi shape=Q N1=10*#meshscale N2=2*#meshscale N3=10*#meshscale N4=2*#meshscale B1=0.5 B2=0.5 B3=0.5 B4=0.5 
/============================================================================================================================
/
/ CLEO detector hollow space. Endpole iron very close to end of solenoid coil
/ Front face of pole tip iron is at 69 inches from IP according to drawing 6056-085.
/ Back face of pole tip is at 105 inches from IP according to drawing 6056-001.
/ Inner radius of pole tip is 18.5/2 inches according to drawings 6052-301 and 6052-303.
/
/ CLEO II and CLEO III
$ const #polendyhi 2.54*105
$ const #detendyhi 69*2.54
$ const #yokexhi 261.
$ const #fixedinrad 95.
/
$if #cleo3 eq 0
/ CLEO II
 $ const #bpradius 2.54*18.5/2.
 $ const #pty1 71*2.54         
 $ const #pty2 #polendyhi         
 $ const #pty3 #polendyhi         
 $ const #ptx1 #bpradius      
 $ const #ptx2 30*2.54      
 $ const #ptx3 #ptx2+(#fixedinrad-#ptx2)/2.      
$else
/ CLEO III  according to drawing 6057-026
 $ const #bpradius 2.54*10.24
 $ const #pty1 94.6*2.54
 $ const #pty2 #pty1         
 $ const #pty3 #polendyhi        
 $ const #ptx1 #bpradius      
 $ const #ptx2 22*2.54     
 $ const #ptx3 28*2.54      
$end if
/
/ ---- Background volume parameters ----
$ const #backheight 2*#polendyhi
$ const #backwidth  2*#yokexhi
$ const #backbase  0.
/
/============================================================================================================================
/ Create yoke
/ The solenoid coil is inside two return yoke layers
/ sandwiching a gap for muon chambers. The yoke layers are 36 cm thick and the gap is 9 cm thick.
/

/ Material 3 - Iron
DRAW +DEFAULTS SHAPE=POLYGON TOLERANCE=5.0E-05 MATERIAL=3 PERM=1000 DENS=0 CONDUCTIVITY=0 PHASE=0 VELOCITY=0 N=0 SYMMETRY=0 XCEN=0 YCEN=0 ANGLE=0 NX=1 DX=0 NY=1 DY=0 ROTATIONS=1 TROTATION=0 MIRROR=NO

/ Set defaults for polygon drawing
POLYGON -DATA +SIZE ESIZE=10*#meshscale  BIAS=0.5 LINESHAPE=STRAIGHT F=NO

/ Two yoke pieces separated by the muon chamber
$const #x1 #yokexhi-36. | $const #y1 -#detendyhi | $const #x2 #yokexhi | $const #y2 #detendyhi
poly +data opt=text x=#x1 y=#y1 | poly +data opt=text x=#x1 y=#y2 | poly +data opt=text x=#x2 y=#y2 | poly +data opt=text x=#x2 y=#y1 | poly opt=close
$const #x1 #yokexhi-81. | $const #y1 -#detendyhi | $const #x2 #yokexhi-45 | $const #y2 #detendyhi
poly +data opt=text x=#x1 y=#y1 | poly +data opt=text x=#x1 y=#y2 | poly +data opt=text x=#x2 y=#y2 | poly +data opt=text x=#x2 y=#y1 | poly opt=close


/ Yoke end pieces connecting the two yoke shells extending down to the pole tip outer radius
$const #x1 #fixedinrad | $const #y1 #detendyhi | $const #x2 #yokexhi | $const #y2 #polendyhi
poly +data opt=text x=#x1 y=#y1 | poly +data opt=text x=#x1 y=#y2 | poly +data opt=text x=#x2 y=#y2 | poly +data opt=text x=#x2 y=#y1 | poly opt=close
yes
$const #x1 #fixedinrad | $const #y1 -#detendyhi | $const #x2 #yokexhi | $const #y2 -#polendyhi
poly +data opt=text x=#x1 y=#y1 | poly +data opt=text x=#x1 y=#y2 | poly +data opt=text x=#x2 y=#y2 | poly +data opt=text x=#x2 y=#y1 | poly opt=close
yes

/ Pole Tips !!!!  Decrease mesh size
POLYGON -DATA +SIZE ESIZE=1*#meshscale  BIAS=0.5 LINESHAPE=STRAIGHT F=NO

poly +data opt=text x=#bpradius y=#detendyhi 
poly +data opt=text x=#ptx1 y=#pty1 
poly +data opt=text x=#ptx2 y=#pty2 
poly +data opt=text x=#ptx3 y=#pty3
poly +data opt=text x=#fixedinrad y=#polendyhi
poly +data opt=text x=#fixedinrad y=#detendyhi 
poly opt=close

poly +data opt=text x=#bpradius y=-#detendyhi 
poly +data opt=text x=#ptx1 y=-#pty1 
poly +data opt=text x=#ptx2 y=-#pty2 
poly +data opt=text x=#ptx3 y=-#pty3
poly +data opt=text x=#fixedinrad y=-#polendyhi
poly +data opt=text x=#fixedinrad y=-#detendyhi 
poly opt=close


/ Air volumes
POLYGON -DATA +SIZE ESIZE=1*#meshscale  BIAS=0.5 LINESHAPE=STRAIGHT F=NO

DRAW +DEFAULTS SHAPE=POLYGON TOLERANCE=5.0E-05 MATERIAL=0 PERM=1 DENS=0 CONDUCTIVITY=0 PHASE=0 VELOCITY=0 N=0 SYMMETRY=0 XCEN=0 YCEN=0 ANGLE=0 NX=1 DX=0 NY=1 DY=0 ROTATIONS=1 TROTATION=0 MIRROR=NO
$const #x1 #bpradius | $const #y1 -#detendyhi | $const #x2 #solradin | $const #y2 #detendyhi
poly +data opt=text x=#x1 y=#y1 F=V V=0 | poly +data opt=text F=NO x=#x1 y=#y2 | poly +data opt=text x=#x2 y=#y2 | poly +data opt=text x=#x2 y=#y1 | poly opt=close
yes
$const #x1 0 | $const #y1 -#polendyhi | $const #x2 #bpradius | $const #y2 #polendyhi
poly +data opt=text x=#x1 y=#y1 F=V V=0 | poly +data opt=text F=NO x=#x1 y=#y2 | poly +data opt=text x=#x2 y=#y2 | poly +data opt=text x=#x2 y=#y1 | poly opt=close
yes
$const #x1 #solradout | $const #y1 -#detendyhi | $const #x2 #yokexhi-81. | $const #y2 #detendyhi
poly +data opt=text x=#x1 y=#y1 F=V V=0 | poly +data opt=text F=NO x=#x1 y=#y2 | poly +data opt=text x=#x2 y=#y2 | poly +data opt=text x=#x2 y=#y1 | poly opt=close
yes
$const #x1 #yokexhi-45. | $const #y1 -#detendyhi | $const #x2 #yokexhi-36. | $const #y2 #detendyhi
poly +data opt=text x=#x1 y=#y1 F=V V=0 | poly +data opt=text F=NO x=#x1 y=#y2 | poly +data opt=text x=#x2 y=#y2 | poly +data opt=text x=#x2 y=#y1 | poly opt=close

/ Fill in pole tips with air to control mesh -- Decrease mesh size
POLYGON -DATA +SIZE ESIZE=1*#meshscale  BIAS=0.5 LINESHAPE=STRAIGHT F=NO

poly +data opt=text x=#ptx1 y=#pty1 
poly +data opt=text x=#ptx2 y=#pty2 
$if #cleo3 eq 1 
 poly +data opt=text x=#ptx3 y=#pty3
 poly +data opt=text x=#ptx2 y=#pty3
$end if
poly +data opt=text x=#ptx1 y=#pty3
poly opt=close

poly +data opt=text x=#ptx1 y=-#pty1 
poly +data opt=text x=#ptx2 y=-#pty2 
$if #cleo3 eq 1 
 poly +data opt=text x=#ptx3 y=-#pty3
 poly +data opt=text x=#ptx2 y=-#pty3
$end if
poly +data opt=text x=#ptx1 y=-#pty3
poly opt=close

/ Define mesh shapes
/CONVERT REG1=1 REG2=* SHAPE=Q

/ Improve mesh
/EDIT REG1=15 REG2=15 OPTION=SIDE NUMBER=1 N=10 BIAS=0.5 CURVATURE=0 F=NO NEIGHBOURS=YES

RECONSTRUCT XMIN=0 XMAX=300 YMIN=-300 YMAX=300 REG1=1 REG2=* MATERIAL=ALL NOT=ANY MESH=NO CLINES=YES BACKGROUND=YES ERASE=YES FILL=REGION NODES=CORNERS LABEL=REGION AXES=YES

$if #interact eq 1
$pause
$end if

/=============================================================================================
/ Create background volume 
/
RECONSTRUCT XMIN=0 XMAX=#backwidth YMIN=-#backheight YMAX=#backheight REG1=1 REG2=* MATERIAL=ALL NOT=ANY MESH=NO CLINES=YES BACKGROUND=YES ERASE=YES FILL=REGION NODES=CORNERS LABEL=REGION AXES=YES

DRAW +DEFAULTS SHAPE=BACKGROUND TOLERANCE=5.0E-05 MATERIAL=0 PERM=1 DENS=0 CONDUCTIVITY=0 PHASE=0 VELOCITY=0 N=0 SYMMETRY=0 XCEN=0 YCEN=0 ANGLE=0 NX=1 DX=0 NY=1 DY=0 ROTATIONS=1 TROTATION=0 MIRROR=NO

POLYGON +DATA OPTION=TEXT X=0   Y=0
POLYGON +DATA OPTION=TEXT X=0   Y=#backheight +SIZE ESIZE=10*#meshscale BIAS=0.1  F=V V=0 LINESHAPE=STRAIGHT 
POLYGON +DATA OPTION=TEXT X=#backwidth Y=0   +SIZE ESIZE=10*#meshscale BIAS=0.01  F=V V=0 LINESHAPE=CURVATURE CURVATURE=1/#backheight
POLYGON +DATA OPTION=TEXT X=0 Y=-#backheight   +SIZE ESIZE=10*#meshscale BIAS=0.99  F=V V=0 LINESHAPE=CURVATURE CURVATURE=1/#backwidth
POLYGON -DATA +SIZE ESIZE=10*#meshscale  BIAS=0.9 LINESHAPE=STRAIGHT F=V V=0
POLYGON OPTION=CLOSE
yes

$if #interact eq 1
$pause
$end if

/Boundary conditions
/EDIT REG1=3 REG2=3 OPTION=SIDE NUMBER=1 F=V V=0 




/=================================================================================================================================
/ Define permeability curve
/=================================================================================================================================
BHDATA MATERIAL=3 TYPE=ISOTROPIC
LOAD 'C:\USER_LOCAL\jac\tentenmerm.bh'
QUIT

/ Generate mesh
MESH +ERRORCHECK TOLERANCE=5.0E-05 -DISPLAY | NO

RECONSTRUCT XMIN=0 XMAX=280 YMIN=-280 YMAX=280 REG1=1 REG2=* MATERIAL=ALL NOT=ANY MESH=YES CLINES=YES BACKGROUND=YES ERASE=YES FILL=REGION NODES=CORNERS LABEL=REGION AXES=REDRAW


$if #interact eq 1
$pause
$end if

/ Define analysis type STATIC
SOLVE TYPE=ST
DATA LINEAR=YES ADIT=0 |QUIT

$os del cleofield_&version&_&cleover&_&current&.op2
/ Write model to file
WRITE FILE=cleofield_&version&_&cleover&_&current&
/Answer yes to overwrite file if already exists
/yes

/ Relax
$OS OPERAANL WRD3=cleofield_&version&_&cleover&_&current& WRD2=ST WRD4=fore
/ Answer yes to close solver window
/yes

READ FILE=C:\USER_LOCAL\jac\2d\cleofield_&version&_&cleover&_&current&.st CASE=1 GEOMETRY=YES


/ Display field
CONTOUR STYLE=LINE AUTOMATIC=YES LINES=50 LABEL=NUMBERS DEFORMED=NO REG1=1 REG2=* MATERIAL=ALL NOT=ANY COMPONENT=POT HOMOGENEITY=NO ERASE=NO

$if #interact eq 1
 $pause
$end if

CONTOUR STYLE=ZONE AUTOMATIC=YES LINES=30 LABEL=NUMBERS DEFORMED=NO REG1=1 REG2=* MATERIAL=ALL NOT=ANY COMPONENT=bmod HOMOGENEITY=NO ERASE=NO

$if #interact eq 1
 $pause
$end if

LINE X1=0 Y1=0 Y2=300 X2=X1 NP=100 CURVATURE=0 COMPONENT=bz HOMOGENEITY=NO ERASE=YES PAUSE=NO AUTOMATIC=YES STYLE=AUTOMATIC COLOUR=AUTOMATIC

$if #integrate eq 1
/ To choose field calculation by integration over currents and magnetizations, use this line
 SET FIELD=MAGNETIC +INTEGRAL FLDTYPE=TOTAL SYMMETRY=AXISYMMETRY SOLUTION=AXI  ELEMENT=LINEAR
$else
 SET FIELD=MAGNETIC  SYMMETRY=AXISYMMETRY SOLUTION=AXI  REFXAXIS=INVE ELEMENT=LINEAR
$end if

$if #table eq 1 

/ Write file of coordinates
$const #nx 180
$const #dx 1
$const #ny 601
$const #dy 0.5
$const #nz 0
$os del points.txt
$open 1 points.txt overwrite
$form 1 string string='           X'
$form 2 string string='           Y'
$form 3 fixed 15 5
$form 4 integer 10
$form 5 string string=' '
$const #np #nx*#ny
$assi 4
$write 1 #np
$assi 1 5 3 5 2 5 3
$do #i 1 #nx
$do #j 1 #ny
$const #x #dx*(#i-1)
$const #y #dy*(#j-1)
$const #z 0.
$write 1 #x #y
$end do
$end do
$close 1

/Read file of coordinates and write field table
$open 1 points.txt read
$read 1 #np
$os del cleofield_&version&_&cleover&_&current&.table
$open 2 cleofield_&version&_&cleover&_&current&.table overwrite
$assi 4 4 4 
$write 2 #nx #ny #nz
$assi 3 3 3 3 3 3
$do #i 1 7
$write 2 #np
$end do
$do #i 1 #np
$read 1 #x #y
point meth=cart xp=#x yp=#y
$write 2 #x #y #z br bz #z
$end do
$close 1 
$close 2

$end if

end
yes