C ------------------------------------------------------------------------ C PROGRAM Magnet C DESIGN MAGNETIC SOLENOIDS THAT GIVE THE SMOOTHEST MAGNETIC C FIELD IN PRESCRIBED REGION. C C OCT.17.1996 C DEVELOPED BY JU LI AT MIT C ------------------------------------------------------------------------ C USE THE FOLLOWING UNITS: C LENGTH UNIT = CM C CURRENT UNIT = AMPERE C MAGNETIC FIELD UNIT = GAUSS C ------------------------------------------------------------------------- PROGRAM Magnet INCLUDE 'm.fh' DIMENSION XX(MAX_OP),F(MAX_OP) PI = dacos(-1.d0) ONE = 1.d0 ZERO = 0.d0 TWO_THIRDS = 2.d0/3.d0 NAME(1) = "X" NAME(2) = "Y" NAME(3) = "Z" NAME(4) = "R" c zero all parameter sets DO 654 IS=1,MAX_SOLENOID DO 654 IV=1,4 MPOLY(IV,IS) = 0 MSINE(IV,IS) = 0 DO 655 I=1,MAX_POLY 655 A(I,IV,IS) = 0.d0 DO 656 I=1,MAX_SINE B(I,IV,IS) = 0.d0 C(I,IV,IS) = 0.d0 D(I,IV,IS) = 0.d0 656 CONTINUE 654 CONTINUE c ===================================================================== c Read instructions from input file "con": READ (*,'(A70)') buf C "Read from configurational file (if 'NULL' then start anew):" READ (*,'(A70)') NAME_READ MAXS = 0 IF (NAME_READ.NE."NULL") THEN c -------------------------------------------------------------- c read configuration from disk OPEN (UNIT=LP_READ, STATUS='UNKNOWN', FORM='FORMATTED', a FILE=NAME_READ) READ (LP_READ,'(A70)') buf C "MAXS = total number of solenoids:" READ (LP_READ,*) MAXS DO 848 IS=1,MAXS READ (LP_READ,'(/,A70,/,A70,/,A70)') buf,buf,buf c "Solenoid No. IS" c '----------------------------------------' c "current density (ampere per cm) of the ISth solenoid:" READ (LP_READ,*) W(IS) DO 843 IV=1,4 READ (LP_READ,'(/,A70)') buf c "MPOLY(IV,IS) = order of polynomials parametrizing X (or Y or Z or R):" READ (LP_READ,*) MPOLY(IV,IS) READ (LP_READ,'(A70)') buf c "order coefficient" DO 845 I=1,MPOLY(IV,IS) 845 READ (LP_READ,*) II, A(II, IV, IS) READ (LP_READ,'(/,A70)') buf c "MSINE(IV,IS) = number of sine functions parametrizing X (or Y or Z or R):" READ (LP_READ,*) MSINE(IV,IS) READ (LP_READ,'(A70)') buf c " bi ki phi_i " DO 842 I=1,MSINE(IV,IS) 842 READ (LP_READ,*) B(I,IV,IS), C(I,IV,IS), D(I,IV,IS) 843 CONTINUE READ (LP_READ,'(A70)') buf c '----------------------------------------' 848 CONTINUE CLOSE(LP_READ) c -------------------------------------------------------------- ENDIF READ (*,'(/,A70)') buf C "Write on configurational file:" READ (*,'(A70)') NAME_WRITE READ (*,'(/,A70)') buf C "How frequently print out the current status (also in "rep"):" READ *,K_OUT c will also write output on this file OPEN (UNIT=LP, STATUS='UNKNOWN', FORM='FORMATTED',FILE='rep') READ (*,'(/,A70)') buf C "How frequently write on configurational file:" READ *,K_WRITE READ (*,'(/,A70)') buf C "(BXMIN BXMAX BYMIN BYMAX BZMIN BZMAX) of the target box:" READ *,BXMIN, BXMAX, BYMIN, BYMAX, BZMIN, BZMAX READ (*,'(/,A70,/,A70,/,A70)') buf,buf,buf C "Number of spatial points in the target box-region that you want C to randomly sample beside default corners(8), face centers(6) and C body center(1):" READ *,NSAMPLE READ (*,'(/,A70)') buf C "Random number generator seed:" READ *,ISEED c initiate random number generator ran1 do i=1,mod(iseed,63346) xxx = ran1(iseed) enddo c generate those targets it = 1 c center tx(1,it) = (bxmin+bxmax)/2. tx(2,it) = (bymin+bymax)/2. tx(3,it) = (bzmin+bzmax)/2. do 603 ix=0,1 do 603 iy=0,1 do 603 iz=0,1 c corners it = it+1 tx(1,it) = bxmin+ix*(bxmax-bxmin) tx(2,it) = bymin+iy*(bymax-bymin) tx(3,it) = bzmin+iz*(bzmax-bzmin) 603 continue do 604 i=it+1,it+6 tx(1,i) = tx(1,1) tx(2,i) = tx(2,1) tx(3,i) = tx(3,1) 604 continue c face centers tx(1,it+1) = tx(1,it+1)-(bxmax-bxmin)/2. tx(1,it+2) = tx(1,it+2)+(bxmax-bxmin)/2. tx(2,it+3) = tx(2,it+3)-(bymax-bymin)/2. tx(2,it+4) = tx(2,it+4)+(bymax-bymin)/2. tx(3,it+5) = tx(3,it+5)-(bzmax-bzmin)/2. tx(3,it+6) = tx(3,it+6)+(bzmax-bzmin)/2. it = it+6 c then the random samplings do 605 i=it+1,it+nsample tx(1,i) = bxmin+ran1(iseed)*(bxmax-bxmin) tx(2,i) = bymin+ran1(iseed)*(bymax-bymin) tx(3,i) = bzmin+ran1(iseed)*(bzmax-bzmin) 605 continue c total number of targets is: nt = it+nsample READ (*,'(/,A70)') buf C "(AXMIN AXMAX AYMIN AYMAX AZMIN AZMAX) of the solenoid box:" READ *,AXMIN, AXMAX, AYMIN, AYMAX, AZMIN, AZMAX MWALL = 6 c Construct geometrical boundary constraints: c wall(1-3) is the normal direction of the wall c wall(4-6) is the coordinate of one point on the wall. do 86 i = 1,mwall do 86 j = 1,6 86 wall(j,i) = 0.d0 wall(1,1) = -1.d0 wall(1,2) = 1.d0 wall(2,3) = -1.d0 wall(2,4) = 1.d0 wall(3,5) = -1.d0 wall(3,6) = 1.d0 wall(4,1) = AXMAX wall(4,2) = AXMIN wall(4,3) = (AXMAX+AXMIN)/2. wall(4,4) = (AXMAX+AXMIN)/2. wall(4,5) = (AXMAX+AXMIN)/2. wall(4,6) = (AXMAX+AXMIN)/2. wall(5,3) = AYMAX wall(5,4) = AYMIN wall(6,5) = AZMAX wall(6,6) = AZMIN READ (*,'(/,A70,/,A70)') buf,buf C "*****************************************" C "Adding how many solenoids:" READ *,n_add DO 9467 IS=MAXS+1,MAXS+N_ADD READ (*,'(/,A70)') buf C "Guess: Current density (ampere per cm) passing the solenoid:" READ *, W(IS) DO IV=1,4 READ (*,'(A70)') buf C "X: order of polynomials number of sines" READ *, MPOLY(IV,IS), MSINE(IV,IS) ENDDO READ (*,'(A70)') buf C "Guess: (X,Y,Z,R) of the first coil:" READ *, A(1,1,IS),A(1,2,IS),A(1,3,IS),A(1,4,IS) READ (*,'(A70)') BUF C "Guess: (X,Y,Z,R) of the last coil:" READ *, XLAST, YLAST, ZLAST, RLAST C do a linear interpolation for the (x,y,z,radius) A(2,1,IS) = XLAST - A(1,1,IS) A(2,2,IS) = YLAST - A(1,2,IS) A(2,3,IS) = ZLAST - A(1,3,IS) A(2,4,IS) = RLAST - A(1,4,IS) c increase the number of solenoids 9467 CONTINUE MAXS = MAXS+n_add READ (*,'(/,A70,/,A70)') buf,buf C "*****************************************" C "Upgrade how many solenoids:" READ *,n_up DO I=1,N_UP READ (*,'(/,A70)') buf C "Solenoid index:" READ *, IS DO IV=1,4 READ (*,'(A70)') buf C "X: order of polynomials number of sines" READ *, MPOLY(IV,IS), MSINE(IV,IS) ENDDO ENDDO c ====================================================================== c assign a (random) reference direction for each solenoid, c as long as it is not parallel with the solenoid axis, c everything will be fine. do is = 1,maxs vx(is) = ran1(iseed) vy(is) = ran1(iseed) vz(is) = ran1(iseed) enddo c Put initial parameter values in XX ig = 0 do 5463 is=1,maxs ig = ig+1 XX(ig) = w(is) do 5463 iv=1,4 do 5454 i=1,MPOLY(IV,IS) ig = ig+1 5454 XX(ig) = A(i,iv,is) do 5463 i=1,MSINE(IV,IS) ig = ig+1 XX(ig) = B(i,iv,is) ig = ig+1 XX(ig) = C(i,iv,is) ig = ig+1 XX(ig) = D(i,iv,is) 5463 CONTINUE c NOP is the total number of optimization parameters NOP = ig call magnetic (1,25.d0,0.d0,0.d0,Bxx,Byy,Bzz) print *,Bxx,Byy,Bzz ITER = 0 CALL FRPRMN(XX,NOP,2.5D-3,ITER,FMIN,VALUE,FORCE) FMIN = DSQRT(FMIN/NT) PRINT *,'THE AVERAGE ERROR CONVERGE TO ',FMIN,' Gauss' PRINT *,'AFTER ',ITER,' ITERATIONS,' PRINT *,'CORRESPONDING TO ',FMIN/5000.*1E6,' ppm.' WRITE (LP,*) 'THE AVERAGE ERROR CONVERGE TO ',FMIN,' Gauss' WRITE (LP,*) 'AFTER ',ITER,' ITERATIONS,' WRITE (LP,*) 'CORRESPONDING TO ',FMIN/5000.*1E6,' ppm.' CALL WRITE_CONFIG CLOSE(LP) STOP END SUBROUTINE WRITE_CONFIG C *************************** C WRITE CONFIGURATIONAL FILE C *************************** INCLUDE 'm.fh' OPEN (UNIT=LP_WRITE, STATUS='UNKNOWN', FORM='FORMATTED', a FILE=NAME_WRITE) WRITE (LP_WRITE,'("Total number of solenoids:")') WRITE (LP_WRITE,'(I2)') MAXS DO 949 IS=1,MAXS WRITE (LP_WRITE,'(/,"Solenoid No.",I2)') IS WRITE (LP_WRITE, a'("------------------------------------------------------")') WRITE (LP_WRITE,'("current density (ampere per cm) of the", I2, a "th solenoid:")') IS WRITE (LP_WRITE,*) W(IS) DO 943 IV=1,4 WRITE (LP_WRITE, '(/, * "order of polynomials parametrizing ",A1)') NAME(IV)(1:1) WRITE (LP_WRITE,'(I2)') MPOLY(IV,IS) WRITE (LP_WRITE,'("order coefficient")') DO 945 I=1,MPOLY(IV,IS) 945 WRITE (LP_WRITE,'(I3,5x,F13.7)') I, A(I, IV, IS) WRITE (LP_WRITE, '(/, * "number of sine functions parametrizing ",A1)') NAME(IV)(1:1) WRITE (LP_WRITE,'(I2)') MSINE(IV,IS) WRITE (LP_WRITE,'(" bi ki phi_i")') DO 942 I=1,MSINE(IV,IS) 942 WRITE (LP_WRITE,'(3(F13.7,2x))')B(I,IV,IS),C(I,IV,IS),D(I,IV,IS) 943 CONTINUE WRITE (LP_WRITE, a'("------------------------------------------------------")') 949 CONTINUE CLOSE(LP_WRITE) RETURN END SUBROUTINE MAGNETIC (IS,X,Y,Z,BX,BY,BZ) C ***************************************************************** C CALCULATE THE MAGNETIC FIELD INDUCED BY SOLENOID No.IS ON C (X,Y,Z) USING ROMBERG INTEGRATION: EPS IS THE DESIRED ACCURACY. C ***************************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (EPS=1.d-6, JMAX=20, JMAXP=JMAX+1, K=5, KM=K-1) DOUBLE PRECISION h(JMAXP),sx(JMAXP),sy(JMAXP),sz(JMAXP) h(1)=1.d0 call coil(IS,0.d0,X,Y,Z,BX0,BY0,BZ0) call coil(IS,1.d0,X,Y,Z,BX1,BY1,BZ1) sx(1) = 0.5*(BX0+BX1) sy(1) = 0.5*(BY0+BY1) sz(1) = 0.5*(BZ0+BZ1) do 12 j=2,JMAX c integration by trapzoid rule with progressively decreasing h h(j)=0.25*h(j-1) it=2**(j-2) tnm=dble(it) del=1.d0/tnm t=0.5*del sumx=0.d0 sumy=0.d0 sumz=0.d0 do jj=1,it call coil(IS,t,X,Y,Z,BXX,BYY,BZZ) sumx=sumx+BXX sumy=sumy+BYY sumz=sumz+BZZ t=t+del enddo sx(j)=0.5*(sx(j-1)+sumx/tnm) sy(j)=0.5*(sy(j-1)+sumy/tnm) sz(j)=0.5*(sz(j-1)+sumz/tnm) if (j.ge.K) then call polint(h(j-KM),sx(j-KM),K,0.d0,BX,dssx) call polint(h(j-KM),sy(j-KM),K,0.d0,BY,dssy) call polint(h(j-KM),sz(j-KM),K,0.d0,BZ,dssz) if (dsqrt(dssx**2+dssy**2+dssz**2).le. a EPS*dsqrt(BX**2+BY**2+BZ**2)) return endif 12 continue pause 'can not converge in Romberg integration, increase EPS' END SUBROUTINE COIL (IS,T,X,Y,Z,BXX,BYY,BZZ) C ************************************************ c calculate the magnetic field induced by a coil c weighted by differential length. C ************************************************ INCLUDE 'a.fh' DOUBLE PRECISION K,K2,K4 CALL axial_point(is,t,xx,yy,zz,radius,cx,cy,cz) c the relative displacement vector of target: rx = x-xx ry = y-yy rz = z-zz r2 = rx*rx+ry*ry+rz*rz r1 = dsqrt(r2) c normalize it erx = rx/r1 ery = ry/r1 erz = rz/r1 c calculate the inner product with the normal direction of the coil cosn = erx*cx+ery*cy+erz*cz ax = erx-cosn*cx ay = ery-cosn*cy az = erz-cosn*cz sinn = ax*erx+ay*ery+az*erz c 'es' stands for "{\bf e}_{\theta}" in my derivation if (sinn.ne.0.d0) then esx = (cosn*erx-cx)/sinn esy = (cosn*ery-cy)/sinn esz = (cosn*erz-cz)/sinn else esx = 0.d0 esy = 0.d0 esz = 0.d0 endif c 'e' stands for "\eta" in my derivation e = r1/radius e2 = e*e de = 1+e2+2*e*sinn k2 = 4*e*sinn/de k = dsqrt(k2) k4 = k2*k2 dkde = k/2.*(1-e2)/e/de c 's' stands for "\theta" dkds = k/2.*cosn*(1+e2)/sinn/de rd1 = rd(ZERO,1-k2,ONE) rd2 = rd(ZERO,ONE,1-k2) rff = rf(ZERO,1-k2,ONE) c 'f' stands for "F(k)" f = rd1*TWO_THIRDS - rff dfdk = (-rd1*TWO_THIRDS+(2-k2)*rd2/3.)/k c 'AA' stands for "A_{\phi}" AA = 0.4*c1*w(is)*F/dsqrt(de) df = 1./k + dfdk/f Br = AA/r1*(cosn/sinn/2.+df*dkds) Bs = -AA/radius*(0.5/e+df*dkde) Bxx = Br*erx+Bs*esx Byy = Br*ery+Bs*esy Bzz = Br*erz+Bs*esz RETURN END DOUBLE PRECISION FUNCTION VALUE(XX) c ********************************************************** c The objective function of conjugate gradient minimization c ********************************************************** INCLUDE 'a.fh' DOUBLE PRECISION XX(MAX_OP) c project XX parameter set to internal variables ig = 0 do 5463 is=1,maxs ig = ig+1 w(is) = XX(ig) do 5463 iv=1,4 do 5454 i=1,MPOLY(IV,IS) ig = ig+1 5454 A(i,iv,is) = XX(ig) do 5463 i=1,MSINE(IV,IS) ig = ig+1 B(i,iv,is) = XX(ig) ig = ig+1 C(i,iv,is) = XX(ig) ig = ig+1 D(i,iv,is) = XX(ig) 5463 CONTINUE c get the magnetic field on all target points do 642 it = 1,nt tbx(it) = 0. tby(it) = 0. tbz(it) = 0. do 642 is = 1,maxs call magnetic(is,tx(1,it),tx(2,it),tx(3,it),bxx,byy,bzz) bx(it,is) = bxx by(it,is) = byy bz(it,is) = bzz tbx(it) = tbx(it)+bxx tby(it) = tby(it)+byy tbz(it) = tbz(it)+bzz 642 continue value = 0.d0 do it = 1,nt c the designed field is -5000 Gauss in z-direction. value = value+tbx(it)**2+tby(it)**2+(tbz(it)+5000)**2 enddo c Evaluate the minimal distance between solenoids, c add repulsions if they touch. Store the derivatives. do 432 isa = 1,maxs do 432 isb = isa+1,maxs dd = distance_solenoids(isa,isb) dist_sw(isa,isb) = dd dist_sw(isb,isa) = dd value = value + repulsion(dd,de) deri_sw(isa,isb) = de deri_sw(isb,isa) = de 432 continue c same for solenoids and walls do is = 1,maxs value = value + penalty_to_walls (is) enddo ITER = ITER+1 if(mod(ITER,K_OUT).eq.0) then print *,'at',ITER,' step, the average error is ', a dsqrt(value),'Gauss' write (LP,5632) ITER, dsqrt(value) 5632 format(I8,' step: error = ',f10.5, 'Gauss') endif RETURN END FUNCTION PENALTY_TO_WALLS (IS) INCLUDE 'm.fh' COMMON /IS/ ISS c pass the common block information: which solenoid ISS = IS c bracket the minima ax = 0. bx = 1. call mnbrak(ax,bx,cx,fa,fb,fc,COIL_TO_WALLS) c golden section search for the minima parametrized by t: c don't need to be very accurate, TOL = 0.001 dist_sw(is) = golden(ax,bx,cx,COIL_TO_WALLS,0.001,tmin) PENALTY_TO_WALLS = REPULSION(dist_sw(is),deri_sw(is)) RETURN END FUNCTION REPULSION (dist,derivative) INCLUDE 'm.fh' PARAMETER (CLOSE_UP=2., PENALTY=5000*5000) if (dist.lt.close_up) then c calculate the imagined repulsion between two objects repulsion = ((1./dist)**12-(1./close_up)**12)*penalty derivative = -12 * (1./dist)**13 * penalty else repulsion = 0. derivative = 0. endif RETURN END FUNCTION COIL_TO_WALLS (T) INCLUDE 'm.fh' COMMON /IS/ ISS coil_to_walls = 1D10 call axial_point (iss,t,x0,y0,z0,radius,cx,cy,cz) do iw = 1,mwall alpha = cx*wall(1,iw)+cy*wall(2,iw)+cz*wall(3,iw) d0 = (x0-wall(4,iw))*wall(1,iw) + (y0-wall(5,iw))*wall(2,iw) a + (z0-wall(6,iw))*wall(3,iw) vd = alpha * radius if ((d0+vd)*(d0-vd).gt.0) then c the coil is on one side of the wall dist = min(dabs(d0+vd),dabs(d0-vd)) coil_to_walls = min(coil_to_walls,dist) else coil_to_walls = 0.d0 return endif enddo RETURN END SUBROUTINE FORCE(XX,F) c ************************************************** c force function of conjugate gradient minimization c ************************************************** INCLUDE 'm.fh' PARAMETER (alpha = 0.001,TOL=1D-5) DOUBLE PRECISION XX(MAX_OP),F(MAX_OP) c Before calling force(), it should have called value(), so c there is no need to re-evaluate XX. c set bound to forces BOUND_W=20.d0 c maximal change in current density = 1000 ampere/cm, etc. BOUND_A=1.d0 BOUND_B=0.5d0 BOUND_C=PI*0.01 BOUND_D=0.1 print *,'evaluating force' ig = 0 do 464 is=1,maxs print *,'solenoid',is if (w(is).gt.tol) then dv = 1. cc = 0. do it=1,nt cc = cc+2*(bx(it,is)/w(is)*tBX(it)+by(it,is)/w(is)*tBY(it)+ a bz(it,is)/w(is)*(tBZ(it)+5000.)) enddo else c if we change the current density dv = BOUND_W*alpha w(is) = w(is)+dv cc = 0. do it=1,nt call magnetic(is,tx(1,it),tx(2,it),tx(3,it),bxx,byy,bzz) cc = cc+2*((bxx-bx(it,is))*tBX(it)+(byy-by(it,is))*tBY(it)+ a (bzz-bz(it,is))*(tBZ(it)+5000.)) enddo w(is) = w(is)-dv endif ig = ig+1 F(ig) = -cc/dv if (dabs(F(ig)).gt.BOUND_W) F(ig) = dsign(BOUND_W, F(ig)) c Do numerical differentiation to get the forces... c if we change one of the following variables do 464 iv=1,4 do 454 i = 1,MPOLY(IV,IS) dv = a(i,iv,is)*alpha if (a(i,iv,is).lt.tol) dv = alpha*BOUND_A a(i,iv,is) = a(i,iv,is) + dv cc = 0. do it=1,nt call magnetic(is,tx(1,it),tx(2,it),tx(3,it),bxx,byy,bzz) cc = cc+2*((bxx-bx(it,is))*tBX(it)+(byy-by(it,is))*tBY(it)+ a (bzz-bz(it,is))*(tBZ(it)+5000.)) enddo a(i,iv,is) = a(i,iv,is) - dv ig = ig+1 F(ig) = -cc/dv if (dabs(F(ig)).gt.BOUND_A) F(ig) = dsign(BOUND_A, F(ig)) 454 continue print *,NAME(iv)(1:1),' polynomials done.' do 463 i = 1,MSINE(IV,IS) dv = b(i,iv,is)*alpha if (b(i,iv,is).lt.tol) dv = alpha*BOUND_B b(i,iv,is) = b(i,iv,is) + dv cc = 0. do it=1,nt call magnetic(is,tx(1,it),tx(2,it),tx(3,it),bxx,byy,bzz) cc = cc+2*((bxx-bx(it,is))*tBX(it)+(byy-by(it,is))*tBY(it)+ a (bzz-bz(it,is))*(tBZ(it)+5000.)) enddo b(i,iv,is) = b(i,iv,is) - dv ig = ig+1 F(ig) = -cc/dv if (dabs(F(ig)).gt.BOUND_B) F(ig) = dsign(BOUND_B, F(ig)) dv = c(i,iv,is)*alpha if (c(i,iv,is).lt.tol) dv = alpha*BOUND_C c(i,iv,is) = c(i,iv,is) + dv cc = 0. do it=1,nt call magnetic(is,tx(1,it),tx(2,it),tx(3,it),bxx,byy,bzz) cc = cc+2*((bxx-bx(it,is))*tBX(it)+(byy-by(it,is))*tBY(it)+ a (bzz-bz(it,is))*(tBZ(it)+5000.)) enddo c(i,iv,is) = c(i,iv,is) - dv ig = ig+1 F(ig) = -cc/dv if (dabs(F(ig)).gt.BOUND_C) F(ig) = dsign(BOUND_C, F(ig)) dv = d(i,iv,is)*alpha if (d(i,iv,is).lt.tol) dv = alpha*BOUND_D d(i,iv,is) = d(i,iv,is) + dv cc = 0. do it=1,nt call magnetic(is,tx(1,it),tx(2,it),tx(3,it),bxx,byy,bzz) cc = cc+2*((bxx-bx(it,is))*tBX(it)+(byy-by(it,is))*tBY(it)+ a (bzz-bz(it,is))*(tBZ(it)+5000.)) enddo d(i,iv,is) = d(i,iv,is) - dv ig = ig+1 F(ig) = -cc/dv if (dabs(F(ig)).gt.BOUND_D) F(ig) = dsign(BOUND_D, F(ig)) 463 continue print *,NAME(iv)(1:1),' sine functions done.' 464 continue RETURN END FUNCTION VALUE_DISTANCE (XX) C *********************************************************** C CALCULATE THE DISTANCE BETWEEN TWO POINTS ON SEPERATE C SOLENOID SURFACES PARAMETRIZED BY (T1,PHI1) AND (T2,PHI2). C *********************************************************** DOUBLE PRECISION XX(4) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /SO/ is1,t1,phi1,is2,t2,phi2,distance,xt1,yt1,zt1, a dx1,dy1,dz1,xt2,yt2,zt2,dx2,dy2,dz2,RX,RY,RZ t1 = xx(1) phi1 = xx(2) t2 = xx(3) phi2 = xx(4) scale = 1.d0 c increase the error if over the boundary if (t1.gt.1.d0) then scale = (t1-1) + scale t1 = 1.d0 elseif (t1.lt.0.d0) then scale = (1-t1) + scale t1 = 0.d0 endif if (t2.gt.1.d0) then scale = (t2-1) + scale t2 = 1.d0 elseif (t2.lt.0.d0) then scale = (1-t2) + scale t2 = 0.d0 endif call surface_point (is1,t1,phi1,xt1,yt1,zt1,dx1,dy1,dz1) call surface_point (is2,t2,phi2,xt2,yt2,zt2,dx2,dy2,dz2) RX = XT1-XT2 RY = YT1-YT2 RZ = ZT1-ZT2 DISTANCE = dsqrt(RX*RX+RY*RY+RZ*RZ) VALUE1 = DISTANCE * scale RETURN END FUNCTION FORCE_DISTANCE (XX,F) C ********************************************************* C CALCULATE HOW MUCH THE DISTANCE WOULD VARY IF WE PERTURB C (T1,PHI1) OR (T2,PHI2). C ********************************************************* IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (VAR=0.001) DIMENSION XX(4),F(4) COMMON /SO/ is1,t1,phi1,is2,t2,phi2,distance,xt1,yt1,zt1, a dx1,dy1,dz1,xt2,yt2,zt2,dx2,dy2,dz2,RX,RY,RZ c if over the boundary then put in negative force if (t1.gt.1.d0) then f(1) = - DISTANCE elseif (t1.lt.0.d0) then f(1) = DISTANCE else c variate t1 to see how much the distance would change dt = dsign(VAR, 0.5-t1) call surface_point (is1,t1+dt,phi1,xt,yt,zt,dx,dy,dz) f(1) = -((xt-xt1)*RX+(yt-yt1)*RY+(zt-zt1)*RZ)/DISTANCE/dt endif c if we variate phi1... can be calculated analytically f(2) = -(RX*DX1+RY*DY1+RZ*DZ1)/DISTANCE c same for t2, phi2 if (t1.gt.1.d0) then f(3) = - DISTANCE elseif (t1.lt.0.d0) then f(3) = DISTANCE else c variate t1 to see how much the distance would change dt = dsign(VAR, 0.5-t2) call surface_point (is2,t2+dt,phi1,xt,yt,zt,dx,dy,dz) f(3) = ((xt-xt2)*RX+(yt-yt2)*RY+(zt-zt2)*RZ)/DISTANCE/dt endif f(4) = (RX*DX2+RY*DY2+RZ*DZ2)/DISTANCE RETURN END FUNCTION DISTANCE_SOLENOIDS (ISA,ISB) C ***************************************************** C CALCULATE THE NEAREST DISTANCE BETWEEN TWO SOLENOIDS C ***************************************************** INCLUDE 'a.fh' COMMON /SO/ is1,t1,phi1,is2,t2,phi2,distance,xt1,yt1,zt1, a dx1,dy1,dz1,xt2,yt2,zt2,dx2,dy2,dz2,RX,RY,RZ DIMENSION XX(4) is1 = isa is2 = isb distance_last = -1. c randomly choose a starting point to ensure global minima 54 xx(1) = ran1(iseed) xx(2) = ran1(iseed)*2*PI xx(3) = ran1(iseed) xx(4) = ran1(iseed)*2*PI call frprmn1 (xx,4,1d-3,iter,fmin) if (dabs(fmin-distance_last).le.1D-3) goto 55 distance_last = fmin goto 54 55 DISTANCE_SOLENOIDS = distance RETURN END SUBROUTINE SURFACE_POINT(IS,T,PHI,XT,YT,ZT,DX,DY,DZ) C *********************************************************** C CALCULATE THE COORDINATES OF A POINT ON SOLENOID SURFACE C #IS, PARAMETRIZED BY (T,PHI). (DX,DY,DZ) IS THE DIFFERENTIAL C DISPLACEMENT IF WE VARY PHI. C *********************************************************** c calculate orgin position of the coil call axial_point (is,t,xx,yy,zz,radius,cx,cy,cz) c use our reference direction to generate an ortho-normal basis call basis(is1,cx,cy,cz,ax,ay,az,bx,by,bz) sinn = dsin(phi) cosn = dcos(phi) XT = XX + ( AX*cosn + BX*sinn ) * radius YT = YY + ( AY*cosn + BY*sinn ) * radius ZT = ZZ + ( AZ*cosn + BZ*sinn ) * radius DX = (-AX*sinn + BX*cosn ) * radius DY = (-AY*sinn + BY*cosn ) * radius DZ = (-AZ*sinn + BZ*cosn ) * radius RETURN END SUBROUTINE AXIAL_POINT (IS,T,XX,YY,ZZ,RADIUS,CX,CY,CZ) C *********************************************************** C CALCULATE THE AXIAL POISTION AND DIRECTION AT PARAMETER T C ON SOLENOID #IS C *********************************************************** INCLUDE 'm.fh' DOUBLE PRECISION TN(MAX_POLY+1) tn(1) = 0.d0 tn(2) = 1.d0 do n=3,MAX_POLY+1 tn(n) = tn(n-1)*t enddo c calculate the position and direction of this coil xx = 0.d0 yy = 0.d0 zz = 0.d0 cx = 0.d0 cy = 0.d0 cz = 0.d0 do i = 1,MPOLY(1,IS) xx = xx+a(i,1,is)*tn(i+1) cx = cx+(i-1)*a(i,1,is)*tn(i) enddo do i = 1,MSINE(1,IS) xx = xx+b(i,1,is)*dsin(c(i,1,is)*(t-d(i,1,is))) cx = cx+b(i,1,is)*c(i,1,is)*dcos(c(i,1,is)*(t-d(i,1,is))) enddo do i = 1,MPOLY(2,IS) yy = yy+a(i,2,is)*tn(i+1) cy = cy+(i-1)*a(i,2,is)*tn(i) enddo do i = 1,MSINE(2,IS) yy = yy+b(i,2,is)*dsin(c(i,2,is)*(t-d(i,2,is))) cy = cy+b(i,2,is)*c(i,2,is)*dcos(c(i,2,is)*(t-d(i,2,is))) enddo do i = 1,MPOLY(3,IS) zz = zz+a(i,3,is)*tn(i+1) cz = cz+(i-1)*a(i,3,is)*tn(i) enddo do i = 1,MSINE(3,IS) zz = zz+b(i,3,is)*dsin(c(i,3,is)*(t-d(i,3,is))) cz = cz+b(i,3,is)*c(i,3,is)*dcos(c(i,3,is)*(t-d(i,3,is))) enddo c calculate the radius of this coil radius = 0.d0 do i = 1,MPOLY(4,IS) radius = radius+a(i,4,is)*tn(i+1) enddo do i = 1,MSINE(4,IS) radius = radius+b(i,4,is)*dsin(c(i,4,is)*(t-d(i,4,is))) enddo c get the normalized direction of the coil c 'c' stand for "{\bf n}" in my derivation c2 = cx*cx+cy*cy+cz*cz c1 = dsqrt(c2) cx = cx/c1 cy = cy/c1 cz = cz/c1 RETURN END SUBROUTINE BASIS (IS,CX,CY,CZ,AX,AY,AZ,BXT,BYT,BZT) C ************************************************* C GENERATE AN MOVING ORTHONORMAL COORDINATE BASIS c USING OUR FIXED RANDOM REFERENCE DIRECTION. C ************************************************* INCLUDE 'm.fh' c take a cross-product V*C AX = VY(IS)*CZ-CY*VZ(IS) AY = CX*VZ(IS)-VX(IS)*CZ AZ = VX(IS)*CY-CX*VY(IS) A1 = dsqrt(AX**2+AY**2+AZ**2) AX = AX/A1 AY = AY/A1 AZ = AZ/A1 c take the cross-product C*A BXT = CY*AZ-AY*CZ BYT = AX*CZ-CX*AZ BZT = CX*AY-AX*CY RETURN END SUBROUTINE polint(xa,ya,n,x,y,dy) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER n,NMAX DOUBLE PRECISION dy,x,y,xa(n),ya(n) PARAMETER (NMAX=10) INTEGER i,m,ns DOUBLE PRECISION den,dif,dift,ho,hp,w,c(NMAX),d(NMAX) ns=1 dif=dabs(x-xa(1)) do 11 i=1,n dift=dabs(x-xa(i)) if (dift.lt.dif) then ns=i dif=dift endif c(i)=ya(i) d(i)=ya(i) 11 continue y=ya(ns) ns=ns-1 do 13 m=1,n-1 do 12 i=1,n-m ho=xa(i)-x hp=xa(i+m)-x w=c(i+1)-d(i) den=ho-hp if(den.eq.0.)pause 'failure in polint' den=w/den d(i)=hp*den c(i)=ho*den 12 continue if (2*ns.lt.n-m)then dy=c(ns+1) else dy=d(ns) ns=ns-1 endif y=y+dy 13 continue return END FUNCTION rd(x,y,z) DOUBLE PRECISION rd,x,y,z,ERRTOL,TINY,BIG,C1,C2,C3,C4,C5,C6 PARAMETER (ERRTOL=.0015d0,TINY=3D-39,BIG=1.7D38,C1=3.d0/14.d0, *C2=1.d0/6.d0,C3=9.d0/22.d0,C4=3.d0/26.d0,C5=.25d0*C3,C6=1.5d0*C4) DOUBLE PRECISION alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac, *sqrtx,sqrty,sqrtz,sum,xt,yt,zt if(min(x,y).lt.0..or.min(x+y,z).lt.TINY.or.max(x,y, *z).gt.BIG)pause 'invalid arguments in rd' xt=x yt=y zt=z sum=0.d0 fac=1.d0 1 continue sqrtx=dsqrt(xt) sqrty=dsqrt(yt) sqrtz=dsqrt(zt) alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz sum=sum+fac/(sqrtz*(zt+alamb)) fac=.25*fac xt=.25*(xt+alamb) yt=.25*(yt+alamb) zt=.25*(zt+alamb) ave=.2*(xt+yt+3.*zt) delx=(ave-xt)/ave dely=(ave-yt)/ave delz=(ave-zt)/ave if(max(dabs(delx),dabs(dely),dabs(delz)).gt.ERRTOL)goto 1 ea=delx*dely eb=delz*delz ec=ea-eb ed=ea-6.*eb ee=ed+ec+ec rd=3.*sum+fac*(1.+ed*(-C1+C5*ed-C6*delz*ee)+delz*(C2*ee+delz*(-C3* *ec+delz*C4*ea)))/(ave*dsqrt(ave)) return END FUNCTION rf(x,y,z) DOUBLE PRECISION rf,x,y,z,ERRTOL,TINY,BIG,THIRD,C1,C2,C3,C4 PARAMETER (ERRTOL=.0025d0,TINY=3d-39,BIG=1.7d38,THIRD=1.d0/3.d0, *C1=1.d0/24.d0,C2=.1d0,C3=3.d0/44.d0,C4=1.d0/14.d0) DOUBLE PRECISION alamb,ave,delx,dely,delz,e2,e3,sqrtx,sqrty, *sqrtz,xt,yt,zt if(min(x,y,z).lt.0..or.min(x+y,x+z,y+z).lt.TINY.or.max(x,y, *z).gt.BIG)pause 'invalid arguments in rf' xt=x yt=y zt=z 1 continue sqrtx=dsqrt(xt) sqrty=dsqrt(yt) sqrtz=dsqrt(zt) alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz xt=.25*(xt+alamb) yt=.25*(yt+alamb) zt=.25*(zt+alamb) ave=THIRD*(xt+yt+zt) delx=(ave-xt)/ave dely=(ave-yt)/ave delz=(ave-zt)/ave if(max(dabs(delx),dabs(dely),dabs(delz)).gt.ERRTOL)goto 1 e2=delx*dely-delz**2 e3=delx*dely*delz rf=(1.+(C1*e2-C2-C3*e3)*e2+C4*e3)/dsqrt(ave) return END FUNCTION ran1(idum) INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV DOUBLE PRECISION ran1,AM,EPS,RNMX PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836, *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS) INTEGER j,k,iv(NTAB),iy SAVE iv,iy DATA iv /NTAB*0/, iy /0/ if (idum.le.0.or.iy.eq.0) then idum=max(-idum,1) do 11 j=NTAB+8,1,-1 k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM j=1+iy/NDIV iy=iv(j) iv(j)=idum ran1=min(AM*iy,RNMX) return END c for the minimization of magnetic field fluctuations INCLUDE 'min.f' c for the minimized distances between solenoids INCLUDE 'min1.f' c for the minimization of distance from solenoid to walls FUNCTION golden(ax,bx,cx,f,tol,xmin) DOUBLE PRECISION golden,ax,bx,cx,tol,xmin,f,R,C EXTERNAL f PARAMETER (R=.61803399,C=1.-R) DOUBLE PRECISION f1,f2,x0,x1,x2,x3 x0=ax x3=cx if(abs(cx-bx).gt.abs(bx-ax))then x1=bx x2=bx+C*(cx-bx) else x2=bx x1=bx-C*(bx-ax) endif f1=f(x1) f2=f(x2) 1 if(abs(x3-x0).gt.tol*(abs(x1)+abs(x2)))then if(f2.lt.f1)then x0=x1 x1=x2 x2=R*x1+C*x3 f1=f2 f2=f(x2) else x3=x2 x2=x1 x1=R*x2+C*x0 f2=f1 f1=f(x1) endif goto 1 endif if(f1.lt.f2)then golden=f1 xmin=x1 else golden=f2 xmin=x2 endif return END