C ------------------------------------------------------------------------ C PROGRAM Magnet (POLYNOMIAL VERSION) C DESIGN MAGNETIC SOLENOIDS THAT GIVE THE SMOOTHEST MAGNETIC C FIELD IN PRESCRIBED REGION. USE DOWNHILL SIMPLEX METHOD. C C NOV.28.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 'mp.fh' DIMENSION XX(MAX_OP) DIMENSION P(MAX_OP+1,MAX_OP),Y(MAX_OP+1),ZZ(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 655 IV=1,4 MPOLY(IV,IS) = 0 DO 655 I=1,MAX_POLY c assign default symmetry to coefficients ON(I,IV,IS) = ((IV.eq.1.or.IV.eq.2.or.IV.eq.4).and.mod(I,2).eq.1) a .or.(IV.eq.3.and.mod(I,2).eq.0) 655 A(I,IV,IS) = 0.d0 654 CONTINUE c ===================================================================== c Read instructions from input file "con": READ (*,'(A70)') buf C "A test session(1) or real job(0):" READ *, IS_TEST READ (*,'(/,A70)') buf C "Read from configurational file (if 'NULL' then start anew):" READ (*,'(A70)') NAME_READ MAXS = 0 IF (NAME_READ(1:4).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 843 I=1,MPOLY(IV,IS) READ (LP_READ,*) II, A(II, IV, IS), ON(II, 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 "Give the initial configuration a perturbation:" READ *, PERTURBATION READ (*,'(/,A70)') buf C "Characteristic current density (ampere per cm) for each solenoid:" READ *,scale_w READ (*,'(/,A70)') buf C "Characteristic length scale of the polynomial coefficients:" READ *,scale_a READ (*,'(/,A70)') buf C "Desired relative accuracy:" READ *,ACCURACY 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(1) tx(1,it) = (bxmin+bxmax)/2. tx(2,it) = (bymin+bymax)/2. tx(3,it) = (bzmin+bzmax)/2. 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(6) 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 do 603 ix=0,1 do 603 iy=0,1 do 603 iz=0,1 c corners(8) 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 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 c the desired magnetic field is -5000 Gauss in z-direction c for each point: do it=1,nt B0(3*it-2) = 0.d0 B0(3*it-1) = 0.d0 B0(3*it) = -5000.d0 enddo READ (*,'(/,A70)') buf C "(AXMIN AXMAX AYMIN AYMAX AZMIN AZMAX) of the solenoid box:" READ *,AXMIN, AXMAX, AYMIN, AYMAX, AZMIN, AZMAX ARMAX = max(AXMAX-AXMIN,AYMAX-AYMIN,AZMAX-AZMIN)/2. 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 "Solenoid #1" DO IV=1,4 READ (*,'(A70)') buf C "X: order of polynomials on/off toggle (high to low) " READ *, MPOLY(IV,IS), FLAG DO I=1,MPOLY(IV,IS) ON(I,IV,IS) = (INT(FLAG).NE.INT(FLAG/10)*10) FLAG = NINT(FLAG/10) ENDDO ENDDO READ (*,'(A70)') buf C "Guess: (X,Y,Z,R) of the first coil:" READ *, XFIRST, YFIRST, ZFIRST, RFIRST 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(1,1,IS) = (XFIRST + XLAST)/2. A(1,2,IS) = (YFIRST + YLAST)/2. A(1,3,IS) = (ZFIRST + ZLAST)/2. A(1,4,IS) = (RFIRST + RLAST)/2. A(2,1,IS) = XLAST - XFIRST A(2,2,IS) = YLAST - YFIRST A(2,3,IS) = ZLAST - ZFIRST A(2,4,IS) = RLAST - RFIRST c assign non-zero initial current density to the solenoid W(IS) = 1.0 9467 CONTINUE c increase the number of solenoids MAXS = MAXS+n_add READ (*,'(/,A70,/,A70)') buf,buf C "*****************************************" C "Upgrade how many solenoids:" READ *,N_UP DO 123 I=1,N_UP READ (*,'(/,A70)') buf C "Solenoid index:" READ *, IS DO 123 IV=1,4 READ (*,'(A70)') buf C "X: order of polynomials on/off toggle (high to low) " READ *, MPOLY(IV,IS), FLAG DO 123 II=1,MPOLY(IV,IS) ON(II,IV,IS) = (INT(FLAG).NE.INT(FLAG/10)*10) FLAG = NINT(FLAG/10) 123 CONTINUE c ====================================================================== c assign a (random) reference direction to 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)-0.5 vy(is) = ran1(iseed)-0.5 vz(is) = ran1(iseed)-0.5 enddo c Put initial parameter values into XX ig = 0 do 5463 is=1,maxs do 5463 iv=1,4 do 5463 i=1,MPOLY(IV,IS) if (on(i,iv,is)) then ig = ig+1 XX(ig) = A(i,iv,is)*(1+perturbation*(ran1(iseed)-0.5)) endif 5463 CONTINUE c NOP is the total number of optimization parameters NOP = ig c ----------------------------------------------------------- print *,'Test session ' print *,' ' call magnetic (1,-149.95d0,0.d0,0.d0,Bxx,Byy,Bzz) print *,Bxx,Byy,Bzz call magnetic (1,-150.d0,0.d0,0.d0,Bxx,Byy,Bzz) print *,Bxx,Byy,Bzz do 640 is = 1,maxs do 640 it = 1,nt call magnetic (is,tx(1,it),tx(2,it),tx(3,it), a sbx(it,is),sby(it,is),sbz(it,is)) 640 continue call calc_magnetic_fluctuation print *,'The initial magnetic fluctuation is', a dsqrt(magnetic_fluctuation),' Gauss.' print *,'frontside B =', tbx(2),tby(2),tbz(2) print *,'box center B =', tbx(1),tby(1),tbz(1) print *,'top B =', tbx(7),tby(7),tbz(7) print *,'backside B =', tbx(3),tby(3),tbz(3) c print *,' ' c print *, distance_to_walls(1) c print *, 'Initial distance between 1st and 2nd solenoid =', c a distance_solenoids(1,2,.false.) print *,' ' print *,'Test complete. ' print *,' ' if (is_test.eq.1) stop c ----------------------------------------------------------- WRITE (*,'(/,"Total number of optimization parameters =",I4)') ig DO 11 IS=1,MAXS WRITE (*,'(/,"Solenoid index #",I4)') IS DO 11 IV=1,4 DO 11 II=1,MPOLY(IV,IS) IF (ON(II,IV,IS)) WRITE (*,'(A1,1X,I2)') name(iv)(1:1), II-1 11 CONTINUE WRITE (*,'(/)') print *,' ' print *,'Constructing vertices for the amoeba.' print *,' ' do 765 ix = 1,nop+1 cc = 0. do ig = 1,nop P(ix,ig) = xx(ig) c construct a random vector zz(ig) = ran1(iseed)-0.5 cc = cc + zz(ig)*zz(ig) enddo c normalize the module cc = dsqrt(cc) do ig=1,nop c add back to the amoeba's vertex if (ix.ne.1) P(ix,ig)=P(ix,ig)+zz(ig)/cc*scale_a zz(ig) = P(ix,ig) enddo y(ix) = value(zz) print *,'Vertex',ix,' error =',dsqrt(y(ix)) 765 continue print *,'Vertex computation complete.' print *,' ' ITER = 0 CALL amoeba(P,y,MAX_OP+1,MAX_OP,NOP,ACCURACY*50,value,niter) print *,' ' c print *,'Amoeba converged. Start conjugate gradient minimization.' c print *,' ' c ig = 0 c do 6463 is=1,maxs c do 6463 iv=1,4 c do 6463 i=1,MPOLY(IV,IS) c if (on(i,iv,is)) then c ig = ig+1 c XX(ig) = A(i,iv,is) c endif c6463 CONTINUE c invoke conjugate gradient minimization procedure CALL FRPRMN(XX,NOP,ACCURACY,ITER,FMIN) FMIN = DSQRT(MAGNETIC_FLUCTUATION) PRINT *,'THE AVERAGE MAGNETIC FLUCTUATION CONVERGE TO ', a FMIN,' Gauss' PRINT *,'AFTER ',NITER,' ITERATIONS,' PRINT *,'CORRESPONDING TO ',FMIN/5000.*1E6,' ppm.' WRITE (LP,*) 'THE AVERAGE MAGNETIC FLUCTUATION CONVERGE TO ' a ,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 'mp.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, a '("current density (ampere per cm) of the",I2,"th solenoid:")') a 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 optimize")') DO 943 I=1,MPOLY(IV,IS) WRITE (LP_WRITE,'(I3,5x,F13.7,5x,L5)') I,A(I,IV,IS),ON(I,IV,IS) c WRITE (LP_WRITE,'(I3,5x,F13.7,5x,A5)') I,A(I,IV,IS),ON(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 #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=25, 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 'mp.fh' DOUBLE PRECISION K,K2,K4 CALL axial_point(is,t,xx,yy,zz,radius,cx,cy,cz,c1) 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) if (r1.ne.0.d0) then c normalize it erx = rx/r1 ery = ry/r1 erz = rz/r1 else c treat as 0+ erx = cx ery = cy erz = cz endif 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 c singular cord -- directly calculate the magentic field c print *,'singular cord' RR = dsqrt(radius*radius + r1*r1) BB = 2 * pi * 0.1 * c1 * w(is) * radius * radius / RR/RR/RR Bxx = BB*cx Byy = BB*cy Bzz = BB*cz return 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,ONE-k2,ONE) rd2 = rd(ZERO,ONE,ONE-k2) rff = rf(ZERO,ONE-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 FUNCTION REPULSION (dist,derivative) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (CLOSE_UP=1., PENALTY=1000, POWER=4) if (dist.lt.close_up) then c set quadratic repulsions between two objects if within close_up repulsion = (1-dist/close_up)**POWER*penalty c derivative smooth derivative = -POWER*(1-dist/close_up)**(POWER-1)/close_up*penalty else repulsion = 0. derivative = 0. endif RETURN END FUNCTION DISTANCE_TO_WALLS (IS) c *********************************************************** c Calculate the nearest distance of a parametrized solenoid c with surrounding walls. If it penetrates then the distance c becomes negative. c *********************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /IS/ ISS EXTERNAL COIL_TO_WALLS c pass the common block information: which solenoid ISS = IS c bracket the minima ax = 0.d0 bx = 1.d0 call mnbrak(ax,bx,cx,fa,fb,fc,COIL_TO_WALLS) c search for the minima parametrized by t, TOL=1D-5 distance_to_walls = brent(ax,bx,cx,COIL_TO_WALLS,1D-5,tmin) RETURN END FUNCTION COIL_TO_WALLS (T) INCLUDE 'mp.fh' COMMON /IS/ ISS DATA ISS /1/ c Make a smooth convex-up interpolation over the boundary. c Because in this minimization we won't use gradient informations, c there is no need to make the 1st-order derivative continuous. if (t.gt.1.) then bonus = (t-1)**2 t = 1. elseif (t.lt.0.) then bonus = t**2 t = 0. else bonus = 0. endif coil_to_walls = 1D10 call axial_point (iss,t,x0,y0,z0,radius,cx,cy,cz,c1) 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 = dsqrt(1-alpha*alpha) * radius c dist is positive inside the wall constraints and negative outside. dist = d0-vd coil_to_walls = min(coil_to_walls,dist) enddo coil_to_walls = coil_to_walls + bonus RETURN END SUBROUTINE CALC_MAGNETIC_FLUCTUATION c **************************************************************** c Solve for the optimal current densities for solenoids under c a specific geometry to get the mean squared magnetic c fluctuation: the input is the specific magnetic fields sbx,y,z. c **************************************************************** INCLUDE 'mp.fh' DIMENSION G(3*MAX_TARGET,MAX_SOLENOID),b(MAX_SOLENOID), a AA(MAX_SOLENOID,MAX_SOLENOID),p(MAX_SOLENOID) nt3 = 3*nt c fill in the G matrix according to sbx,y,z do 45 j=1,maxs do 45 i=1,nt G(3*i-2,j)=sbx(i,j)/w(j) G(3*i-1,j)=sby(i,j)/w(j) G(3*i ,j)=sbz(i,j)/w(j) 45 continue c Construct AA=G'*G+nt/(scale_w)^2/maxs: c only the upper triangle is needed. do 46 i=1,maxs do 46 j=i,maxs AA(i,j) = 0.d0 do 47 k=1,nt3 47 AA(i,j)=AA(i,j)+G(k,i)*G(k,j) if (i.eq.j) AA(i,j)=AA(i,j)+nt/scale_w/scale_w/maxs 46 continue c construct b=G'*B0: do 48 i=1,maxs b(i)=0.d0 do 49 k=1,nt3 49 b(i)=b(i)+G(k,i)*B0(k) 48 continue c Cholesky solver of linear equations for symmetric and c positive definite matrix: AA*w=b call choldc(AA,maxs,MAX_SOLENOID,p) call cholsl(AA,maxs,MAX_SOLENOID,p,b,w) c refresh sbx,y,z and get magnetic field at each target point do 75 i=1,nt tbx(i) = 0. tby(i) = 0. tbz(i) = 0. do 75 j=1,maxs sbx(i,j)=G(3*i-2,j)*w(j) sby(i,j)=G(3*i-1,j)*w(j) sbz(i,j)=G(3*i,j)*w(j) tbx(i)=tbx(i)+sbx(i,j) tby(i)=tby(i)+sby(i,j) tbz(i)=tbz(i)+sbz(i,j) 75 continue c get the mean squared magnetic fluctuation magnetic_fluctuation = 0.d0 do 77 i=1,nt magnetic_fluctuation = magnetic_fluctuation+ a(tbx(i)-b0(3*i-2))**2+(tby(i)-b0(3*i-1))**2+(tbz(i)-b0(3*i))**2 77 continue magnetic_fluctuation = magnetic_fluctuation/nt c In order to prevent too big current densities, add penalty c function according to its characteristic scale. current_compensation = 0.d0 do 79 j=1,maxs 79 current_compensation = current_compensation+(w(j)/scale_w)**2 current_compensation = current_compensation/maxs total_fluctuation = magnetic_fluctuation + current_compensation RETURN END DOUBLE PRECISION FUNCTION VALUE(XX) c *************************************** c The objective function of minimization c *************************************** INCLUDE 'mp.fh' DOUBLE PRECISION XX(MAX_OP) c project XX parameters to internal variables ig = 0 do 5463 is=1,maxs do 5463 iv=1,4 do 5463 i=1,MPOLY(IV,IS) if (on(i,iv,is)) then ig = ig+1 A(i,iv,is) = XX(ig) endif 5463 CONTINUE value = 0. touch = .false. c ------------------------------------------------- c calculate the nearest distance between solenoids, c add repulsion if they touch, store the derivatives. do 432 isa = 1,maxs do 432 isb = isa+1,maxs dd = distance_solenoids(isa,isb,.FALSE.) c start_from_last = false: search for global minima. dist_ss(isa,isb) = dd dist_ss(isb,isa) = dd if (dd.le.0.d0) touch = .true. value = value + repulsion(dd,de) deri_ss(isa,isb) = de deri_ss(isb,isa) = de 432 continue c same for solenoids and walls do is = 1,maxs dd = distance_to_walls (is) if (dd.le.0.d0) touch = .true. dist_sw(is) = dd value = value + repulsion(dd,deri_sw(is)) enddo c ------------------------------------------------- if (touch) then print *,' TOUCH: value = ', value*2 c there is no need to calculate the magnetics -- c we construct a plateau with scale factor 2. value = value*2 return endif c calculate the specific magnetic fields sbx,sby,sbz. do 642 it = 1,nt do 642 is = 1,maxs if (w(is).eq.0.) w(is) = 1.d0 call magnetic (is,tx(1,it),tx(2,it),tx(3,it), a sbx(it,is),sby(it,is),sbz(it,is)) 642 continue c calculate the minimum magnetic fluctuation at this geometry call calc_magnetic_fluctuation value = value + total_fluctuation iter = iter+1 if(mod(iter,k_out).eq.0) then print *,'at',iter,' step, the average error is ', a dsqrt(value) print *,'in which',dsqrt(magnetic_fluctuation), a 'Gauss is due to magnetic fluctuation.' write (LP,5632) ITER, dsqrt(value),dsqrt(magnetic_fluctuation) 5632 format(I8,' step: error = ',f10.5,'; magnetic fluctuation = ', a f10.5,' Gauss.') endif RETURN END SUBROUTINE FORCE(XX,F) c ************************************** c force field of the objective function c ************************************** INCLUDE 'mp.fh' PARAMETER (alpha=0.001) DOUBLE PRECISION XX(MAX_OP),F(MAX_OP) c Before calling force(), it should have called value(), so c there is no need to evaluate XX. print *,'evaluating force' old_total_fluctuation = total_fluctuation ig = 0 do 464 is=1,maxs print *,'solenoid no.',is do 464 iv=1,4 do 454 i = 1,MPOLY(IV,IS) if (dabs(a(i,iv,is)).lt.SCALE_A) then dv = alpha*SCALE_A else dv = alpha*a(i,iv,is) endif a(i,iv,is) = a(i,iv,is) + dv cc = 0. c ---------------------------------------------------------- c Evaluate the forces due to artificial repulsion: c solenoid-solenoid do 651 isb = 1,maxs if (is.eq.isb) goto 651 dist = distance_solenoids(is,isb,.TRUE.) c start_from_last = true: start from the previous minima c position between the two solenoids. cc = cc+deri_ss(is,isb)*(dist-dist_ss(is,isb)) 651 continue c solenoid-wall cc = cc+deri_sw(is)*(distance_to_walls(is)-dist_sw(is)) c ---------------------------------------------------------- c assume first that it's on a plateau cc = cc*2 if (.not.touch) then c ------------------------------------------------------------ if (w(is).eq.0) w(is) = 1.d0 do it=1,nt c re-evaluate the specific magnetic fields of this solenoid: call magnetic (is,tx(1,it),tx(2,it),tx(3,it), a sbx(it,is),sby(it,is),sbz(it,is)) enddo call calc_magnetic_fluctuation cc = cc/2 + total_fluctuation - old_total_fluctuation c ------------------------------------------------------------- endif a(i,iv,is) = a(i,iv,is) - dv ig = ig+1 F(ig) = -cc/dv 454 continue print *,NAME(iv)(1:1),' polynomials 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 *********************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION XX(4) PARAMETER (TINY=5D-5) COMMON /SO/ is1,is2,t1,phi1,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) bonus = 0.d0 c Increase the distance if over the boundary, this time we c need to construct a local minima with smooth 1st-derivative. if (t1.gt.1.d0) then bonus = (t1-1-tiny)**2 - tiny**2 + bonus t1 = 1.d0 elseif (t1.lt.0.d0) then bonus = (t1+tiny)**2 - tiny**2 + bonus t1 = 0.d0 endif if (t2.gt.1.d0) then bonus = (t2-1-tiny)**2 - tiny**2 + bonus t2 = 1.d0 elseif (t2.lt.0.d0) then bonus = (t2+tiny)**2 - tiny**2 + bonus 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) c now it is always positive... value_distance = distance + bonus RETURN END SUBROUTINE FORCE_DISTANCE (XX,F) C ********************************************************* C CALCULATE HOW MUCH THE DISTANCE WOULD VARY IF WE PERTURB C (T1,PHI1) OR (T2,PHI2). C ********************************************************* INCLUDE 'mp.fh' PARAMETER (VAR=0.001) DIMENSION XX(4),F(4) PARAMETER (TINY=5D-5) COMMON /SO/ is1,is2,t1,phi1,t2,phi2,distance,xt1,yt1,zt1, a dx1,dy1,dz1,xt2,yt2,zt2,dx2,dy2,dz2,RX,RY,RZ data distance /100./ t1 = xx(1) phi1 = xx(2) t2 = xx(3) phi2 = xx(4) c if over the boundary then put in negative force if (t1.gt.1.d0) then f(1) = -2 * (t1-1-tiny) elseif (t1.lt.0.d0) then f(1) = -2 * (t1+tiny) else c numerically 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 (t2.gt.1.d0) then f(3) = - 2 * (t2-1-tiny) elseif (t2.lt.0.d0) then f(3) = - 2 * (t2+tiny) else dt = dsign(var, 0.5-t2) call surface_point (is2,t2+dt,phi2,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,START_FROM_LAST) C ***************************************************** C CALCULATE THE NEAREST DISTANCE BETWEEN TWO SOLENOIDS C ***************************************************** INCLUDE 'mp.fh' PARAMETER (TINY=5D-5) LOGICAL START_FROM_LAST COMMON /SO/ is1,is2,t1,phi1,t2,phi2,distance,xt1,yt1,zt1, a dx1,dy1,dz1,xt2,yt2,zt2,dx2,dy2,dz2,RX,RY,RZ DIMENSION XX(4),XX_LAST(4,MAX_SOLENOID,MAX_SOLENOID) SAVE XX_LAST c ----------------------------------------------------------------- c If START_FROM_LAST is true, then we are doing a force evaluation, c then necessarily we should start from the previous minimized c configuration. c c If START_FROM_LAST is false, then we need to save the configuration c at the end of the minimization. Furthermore, to ensure global c convergence, we need to start with a random point and re-check. c ------------------------------------------------------------------ is1 = isa is2 = isb distance_last = 1D10 54 if (.not.START_FROM_LAST) then c Randomly pick a starting point, xx(1) = ran1(iseed) xx(2) = ran1(iseed)*2*PI xx(3) = ran1(iseed) xx(4) = ran1(iseed)*2*PI else c start with the previous minimized configuation. xx(1) = xx_last(1,isa,isb) xx(2) = xx_last(2,isa,isb) xx(3) = xx_last(3,isa,isb) xx(4) = xx_last(4,isa,isb) endif iter1 = 0 call frprmn1 (xx,4,TINY*1D-3,iter1,fmin) if (dabs(distance-distance_last).le.TINY.or.distance.le.TINY) a goto 55 distance_last = min(distance, distance_last) c re-check to see if indeed is the global minima. if (.not.START_FROM_LAST) goto 54 55 if (distance.le.TINY) then c the solenoids are touching, estimate how much in depth c ==> negative value to distance_solenoids call axial_point (is1,xx(1),x01,y01,z01,r1,cx1,cy1,cz1,c11) c calculate normal direction of the solenoid surface. rx1 = xt1 - x01 ry1 = yt1 - y01 rz1 = zt1 - z01 call axial_point (is2,xx(3),x02,y02,z02,r2,cx2,cy2,cz2,c21) rx2 = xt2 - x02 ry2 = yt2 - y02 rz2 = zt2 - z02 distance_solenoids = min(r1,r2) * a (((rx1*rx2+ry1*ry2+rz1*rz2)/r1/r2)**2-1) else distance_solenoids = distance endif if (.not.START_FROM_LAST) then c save the current minima position for future force evaluation xx_last(1,isa,isb) = xx(1) xx_last(2,isa,isb) = xx(2) xx_last(3,isa,isb) = xx(3) xx_last(4,isa,isb) = xx(4) xx_last(1,isb,isa) = xx(1) xx_last(2,isb,isa) = xx(2) xx_last(3,isb,isa) = xx(3) xx_last(4,isb,isa) = xx(4) endif c print *, 'the true distance is ',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 *********************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Z) c calculate orgin position of the coil call axial_point (is,t,xx,yy,zz,radius,cx,cy,cz,c1) c use our reference direction to generate an orthonormal basis call basis(is,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,C1) C *********************************************************** C CALCULATE THE AXIAL POISTION AND DIRECTION AT PARAMETER T C ON SOLENOID #IS C *********************************************************** INCLUDE 'mp.fh' DOUBLE PRECISION TN(MAX_POLY+1) c shift the center to t=0.5: better symmetry. tn(1) = 0.d0 tn(2) = 1.d0 do n=3,MAX_POLY+1 tn(n) = tn(n-1)*(t-0.5) 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,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,MPOLY(3,IS) zz = zz+a(i,3,is)*tn(i+1) cz = cz+(i-1)*a(i,3,is)*tn(i) 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 radius = dabs(radius) 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,BX,BY,BZ) include 'mp.fh' C ************************************************* C GENERATE AN MOBILE ORTHONORMAL COORDINATE BASIS c USING OUR FIXED RANDOM REFERENCE DIRECTION. C ************************************************* 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 BX = CY*AZ-AY*CZ BY = AX*CZ-CX*AZ BZ = 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) then pause 'invalid arguments in rd' endif 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 c INCLUDE 'min.f' SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,func) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION ax,bx,cx,fa,fb,fc,func,GOLD,GLIMIT,TINY EXTERNAL func PARAMETER (GOLD=1.618034, GLIMIT=10.d0, TINY=1.d-20) DOUBLE PRECISION dum,fu,q,r,u,ulim fa=func(ax) fb=func(bx) if(fb.gt.fa)then dum=ax ax=bx bx=dum dum=fb fb=fa fa=dum endif cx=bx+GOLD*(bx-ax) fc=func(cx) 1 if(fb.ge.fc)then r=(bx-ax)*(fb-fc) q=(bx-cx)*(fb-fa) u=bx-((bx-cx)*q-(bx-ax)*r)/(2.d0*sign(max(dabs(q-r),TINY),q-r)) ulim=bx+GLIMIT*(cx-bx) if((bx-u)*(u-cx).gt.0.d0)then fu=func(u) if(fu.lt.fc)then ax=bx fa=fb bx=u fb=fu return else if(fu.gt.fb)then cx=u fc=fu return endif u=cx+GOLD*(cx-bx) fu=func(u) else if((cx-u)*(u-ulim).gt.0.d0)then fu=func(u) if(fu.lt.fc)then bx=cx cx=u u=cx+GOLD*(cx-bx) fb=fc fc=fu fu=func(u) endif else if((u-ulim)*(ulim-cx).ge.0.d0)then u=ulim fu=func(u) else u=cx+GOLD*(cx-bx) fu=func(u) endif ax=bx bx=cx cx=u fa=fb fb=fc fc=fu goto 1 endif return END DOUBLE PRECISION FUNCTION brent(ax,bx,cx,f,tol,xmin) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER ITMAX DOUBLE PRECISION ax,bx,cx,tol,xmin,f,CGOLD,ZEPS EXTERNAL f PARAMETER (ITMAX=1000,CGOLD=.3819660,ZEPS=1.0d-10) INTEGER iterr DOUBLE PRECISION a,b,d,e,etemp,fu,fv,fw,fx,p,q,r, a tol1,tol2,u,v,w,x,xm a=min(ax,cx) b=max(ax,cx) v=bx w=v x=v e=0.d0 fx=f(x) fv=fx fw=fx do 11 iterr=1,ITMAX xm=0.5*(a+b) tol1=tol*dabs(x)+ZEPS tol2=2.d0*tol1 if(dabs(x-xm).le.(tol2-.5*(b-a))) goto 3 if(dabs(e).gt.tol1) then r=(x-w)*(fx-fv) q=(x-v)*(fx-fw) p=(x-v)*q-(x-w)*r q=2.*(q-r) if(q.gt.0.) p=-p q=dabs(q) etemp=e e=d if(dabs(p).ge.dabs(.5*q*etemp).or.p.le.q*(a-x).or.p.ge.q*(b-x)) *goto 1 d=p/q u=x+d if(u-a.lt.tol2 .or. b-u.lt.tol2) d=sign(tol1,xm-x) goto 2 endif 1 if(x.ge.xm) then e=a-x else e=b-x endif d=CGOLD*e 2 if(dabs(d).ge.tol1) then u=x+d else u=x+sign(tol1,d) endif fu=f(u) if(fu.le.fx) then if(u.ge.x) then a=x else b=x endif v=w fv=fw w=x fw=fx x=u fx=fu else if(u.lt.x) then a=u else b=u endif if(fu.le.fw .or. w.eq.x) then v=w fv=fw w=u fw=fu else if(fu.le.fv .or. v.eq.x .or. v.eq.w) then v=u fv=fu endif endif 11 continue pause 'brent exceed maximum iterations' 3 xmin=x brent=fx return END DOUBLE PRECISION FUNCTION f1dim(x) include 'mp.fh' DOUBLE PRECISION VALUE,x CU USES VALUE INTEGER j,ncom DOUBLE PRECISION pcom(MAX_OP),xicom(MAX_OP),xt(MAX_OP) COMMON /f1com/ pcom,xicom,ncom do 11 j=1,ncom xt(j)=pcom(j)+x*xicom(j) 11 continue f1dim=VALUE(xt) return END SUBROUTINE linmin(p,xi,n,fret) include 'mp.fh' INTEGER n DOUBLE PRECISION fret,p(n),xi(n),TOL PARAMETER (TOL=1.d-2) CU USES brent,f1dim,mnbrak INTEGER j,ncom DOUBLE PRECISION ax,bx,fa,fb,fx,xmin,xx,pcom(MAX_OP), a xicom(MAX_OP),brent COMMON /f1com/ pcom,xicom,ncom EXTERNAL f1dim ncom=n do 11 j=1,n pcom(j)=p(j) xicom(j)=xi(j) 11 continue ax=0. xx=1. call mnbrak(ax,xx,bx,fa,fx,fb,f1dim) fret=brent(ax,xx,bx,f1dim,TOL,xmin) do 12 j=1,n xi(j)=xmin*xi(j) p(j)=p(j)+xi(j) 12 continue return END SUBROUTINE frprmn(p,n,ftol,iterr,fret) include 'mp.fh' INTEGER iterr,n,ITMAX DOUBLE PRECISION fret,ftol,p(n),EPS,VALUE EXTERNAL VALUE PARAMETER (ITMAX=100000,EPS=1.d-7) CU USES FORCE,VALUE,linmin INTEGER its,j DOUBLE PRECISION dgg,fp,gam,gg,g(MAX_OP),h(MAX_OP),xi(MAX_OP) fp=VALUE(p) call FORCE(p,xi) do 11 j=1,n g(j)=-xi(j) h(j)=g(j) xi(j)=h(j) 11 continue do 14 its=1,ITMAX iterr=its call linmin(p,xi,n,fret) if(2.*dabs(fret-fp).le.ftol*(dabs(fret)+dabs(fp)+EPS)) return fp=VALUE(p) call FORCE(p,xi) gg=0. dgg=0. do 12 j=1,n gg=gg+g(j)**2 C dgg=dgg+xi(j)**2 dgg=dgg+(xi(j)+g(j))*xi(j) 12 continue if(gg.eq.0.)return gam=dgg/gg do 13 j=1,n g(j)=-xi(j) h(j)=g(j)+gam*h(j) xi(j)=h(j) 13 continue 14 continue pause 'frprmn maximum iterations exceeded' return END FUNCTION amotry(p,y,psum,mp,np,ndim,funk,ihi,fac) INTEGER ihi,mp,ndim,np,NMAX DOUBLE PRECISION amotry,fac,p(mp,np),psum(np),y(mp),funk PARAMETER (NMAX=1000) EXTERNAL funk CU USES funk INTEGER j DOUBLE PRECISION fac1,fac2,ytry,ptry(NMAX) fac1=(1.-fac)/ndim fac2=fac1-fac do 11 j=1,ndim ptry(j)=psum(j)*fac1-p(ihi,j)*fac2 11 continue ytry=funk(ptry) if (ytry.lt.y(ihi)) then y(ihi)=ytry do 12 j=1,ndim psum(j)=psum(j)-p(ihi,j)+ptry(j) p(ihi,j)=ptry(j) 12 continue endif amotry=ytry return END SUBROUTINE amoeba(p,y,mp,np,ndim,ftol,funk,iter) INTEGER iter,mp,ndim,np,NMAX,ITMAX DOUBLE PRECISION ftol,p(mp,np),y(mp),funk PARAMETER (NMAX=1000,ITMAX=5000) EXTERNAL funk CU USES amotry,funk INTEGER i,ihi,ilo,inhi,j,m,n DOUBLE PRECISION rtol,sum,swap,ysave,ytry,psum(NMAX),amotry iter=0 1 do 12 n=1,ndim sum=0. do 11 m=1,ndim+1 sum=sum+p(m,n) 11 continue psum(n)=sum 12 continue 2 ilo=1 if (y(1).gt.y(2)) then ihi=1 inhi=2 else ihi=2 inhi=1 endif do 13 i=1,ndim+1 if(y(i).le.y(ilo)) ilo=i if(y(i).gt.y(ihi)) then inhi=ihi ihi=i else if(y(i).gt.y(inhi)) then if(i.ne.ihi) inhi=i endif 13 continue rtol=2.*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo))) if (rtol.lt.ftol) then swap=y(1) y(1)=y(ilo) y(ilo)=swap do 14 n=1,ndim swap=p(1,n) p(1,n)=p(ilo,n) p(ilo,n)=swap 14 continue return endif if (iter.ge.ITMAX) pause 'ITMAX exceeded in amoeba' iter=iter+2 ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,-1.d0) if (ytry.le.y(ilo)) then ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,2.d0) else if (ytry.ge.y(inhi)) then ysave=y(ihi) ytry=amotry(p,y,psum,mp,np,ndim,funk,ihi,0.d5) if (ytry.ge.ysave) then do 16 i=1,ndim+1 if(i.ne.ilo)then do 15 j=1,ndim psum(j)=0.5*(p(i,j)+p(ilo,j)) p(i,j)=psum(j) 15 continue y(i)=funk(psum) endif 16 continue iter=iter+ndim goto 1 endif else iter=iter-1 endif goto 2 END c for the minimized distances between solenoids c INCLUDE 'min1.f' SUBROUTINE mnbrak1 (ax,bx,cx,fa,fb,fc,func) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION ax,bx,cx,fa,fb,fc,func,GOLD,GLIMIT,TINY EXTERNAL func PARAMETER (GOLD=1.618034, GLIMIT=10.d0, TINY=1.d-20) DOUBLE PRECISION dum,fu,q,r,u,ulim fa=func(ax) fb=func(bx) if(fb.gt.fa)then dum=ax ax=bx bx=dum dum=fb fb=fa fa=dum endif cx=bx+GOLD*(bx-ax) fc=func(cx) 1 if(fb.ge.fc)then r=(bx-ax)*(fb-fc) q=(bx-cx)*(fb-fa) u=bx-((bx-cx)*q-(bx-ax)*r)/(2.d0*sign(max(dabs(q-r),TINY),q-r)) ulim=bx+GLIMIT*(cx-bx) if((bx-u)*(u-cx).gt.0.d0)then fu=func(u) if(fu.lt.fc)then ax=bx fa=fb bx=u fb=fu return else if(fu.gt.fb)then cx=u fc=fu return endif u=cx+GOLD*(cx-bx) fu=func(u) else if((cx-u)*(u-ulim).gt.0.d0)then fu=func(u) if(fu.lt.fc)then bx=cx cx=u u=cx+GOLD*(cx-bx) fb=fc fc=fu fu=func(u) endif else if((u-ulim)*(ulim-cx).ge.0.d0)then u=ulim fu=func(u) else u=cx+GOLD*(cx-bx) fu=func(u) endif ax=bx bx=cx cx=u fa=fb fb=fc fc=fu goto 1 endif return END DOUBLE PRECISION FUNCTION brent1 (ax,bx,cx,f,tol,xmin) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER ITMAX DOUBLE PRECISION ax,bx,cx,tol,xmin,f,CGOLD,ZEPS EXTERNAL f PARAMETER (ITMAX=1000,CGOLD=.3819660,ZEPS=1.0d-10) INTEGER iter DOUBLE PRECISION a,b,d,e,etemp,fu,fv,fw,fx,p,q,r, a tol1,tol2,u,v,w,x,xm a=min(ax,cx) b=max(ax,cx) v=bx w=v x=v e=0.d0 fx=f(x) fv=fx fw=fx do 11 iter=1,ITMAX xm=0.5*(a+b) tol1=tol*dabs(x)+ZEPS tol2=2.d0*tol1 if(dabs(x-xm).le.(tol2-.5*(b-a))) goto 3 if(dabs(e).gt.tol1) then r=(x-w)*(fx-fv) q=(x-v)*(fx-fw) p=(x-v)*q-(x-w)*r q=2.*(q-r) if(q.gt.0.) p=-p q=dabs(q) etemp=e e=d if(dabs(p).ge.dabs(.5*q*etemp).or.p.le.q*(a-x).or.p.ge.q*(b-x)) *goto 1 d=p/q u=x+d if(u-a.lt.tol2 .or. b-u.lt.tol2) d=sign(tol1,xm-x) goto 2 endif 1 if(x.ge.xm) then e=a-x else e=b-x endif d=CGOLD*e 2 if(dabs(d).ge.tol1) then u=x+d else u=x+sign(tol1,d) endif fu=f(u) if(fu.le.fx) then if(u.ge.x) then a=x else b=x endif v=w fv=fw w=x fw=fx x=u fx=fu else if(u.lt.x) then a=u else b=u endif if(fu.le.fw .or. w.eq.x) then v=w fv=fw w=u fw=fu else if(fu.le.fv .or. v.eq.x .or. v.eq.w) then v=u fv=fu endif endif 11 continue pause 'brent1 exceed maximum iterations' 3 xmin=x brent1=fx return END DOUBLE PRECISION FUNCTION f1dim1 (x) DOUBLE PRECISION VALUE_DISTANCE,x EXTERNAL VALUE_DISTANCE CU USES VALUE_DISTANCE INTEGER j,ncom1 DOUBLE PRECISION pcom1(4),xicom1(4),xt(4) COMMON /f1com1/ pcom1,xicom1,ncom1 do 11 j=1,ncom1 xt(j)=pcom1(j)+x*xicom1(j) 11 continue f1dim1 = VALUE_DISTANCE(xt) return END SUBROUTINE linmin1 (p,xi,n,fret) INTEGER n DOUBLE PRECISION fret,p(4),xi(4),TOL PARAMETER (TOL=1.d-7) CU USES brent1,f1dim1,mnbrak1 INTEGER j,ncom1 DOUBLE PRECISION ax,bx,fa,fb,fx,xmin,xx,pcom1(4), a xicom1(4),brent1 COMMON /f1com1/ pcom1,xicom1,ncom1 EXTERNAL f1dim1 ncom1=n do 11 j=1,n pcom1(j)=p(j) xicom1(j)=xi(j) 11 continue ax=0. xx=1. call mnbrak1(ax,xx,bx,fa,fx,fb,f1dim1) fret=brent1(ax,xx,bx,f1dim1,TOL,xmin) do 12 j=1,n xi(j)=xmin*xi(j) p(j)=p(j)+xi(j) 12 continue return END SUBROUTINE frprmn1 (p,n,ftol,iter,fret) INTEGER iter,n,ITMAX DOUBLE PRECISION fret,ftol,p(4),EPS,VALUE_DISTANCE EXTERNAL VALUE_DISTANCE PARAMETER (ITMAX=10000,EPS=1.d-7) CU USES FORCE_DISTANCE,VALUE_DISTANCE,linmin1 INTEGER its,j DOUBLE PRECISION dgg,fp,gam,gg,g(4),h(4),xi(4) fp=VALUE_DISTANCE(p) call FORCE_DISTANCE(p,xi) do 11 j=1,n g(j)=-xi(j) h(j)=g(j) xi(j)=h(j) 11 continue do 14 its=1,ITMAX iter=its call linmin1(p,xi,n,fret) if(2.*dabs(fret-fp).le.ftol*(dabs(fret)+dabs(fp)+EPS)) return fp=VALUE_DISTANCE(p) call FORCE_DISTANCE(p,xi) gg=0. dgg=0. do 12 j=1,n gg=gg+g(j)**2 C dgg=dgg+xi(j)**2 dgg=dgg+(xi(j)+g(j))*xi(j) 12 continue if(gg.eq.0.)return gam=dgg/gg do 13 j=1,n g(j)=-xi(j) h(j)=g(j)+gam*h(j) xi(j)=h(j) 13 continue 14 continue pause 'frprmn1 maximum iterations exceeded' return END SUBROUTINE choldc(a,n,np,p) INTEGER n,np DOUBLE PRECISION a(np,np),p(n) INTEGER i,j,k DOUBLE PRECISION sum do 13 i=1,n do 12 j=i,n sum=a(i,j) do 11 k=i-1,1,-1 sum=sum-a(i,k)*a(j,k) 11 continue if(i.eq.j)then if(sum.le.0.)pause 'choldc failed' p(i)=sqrt(sum) else a(j,i)=sum/p(i) endif 12 continue 13 continue return END SUBROUTINE cholsl(a,n,np,p,b,x) INTEGER n,np DOUBLE PRECISION a(np,np),b(n),p(n),x(n) INTEGER i,k DOUBLE PRECISION sum do 12 i=1,n sum=b(i) do 11 k=i-1,1,-1 sum=sum-a(i,k)*x(k) 11 continue x(i)=sum/p(i) 12 continue do 14 i=n,1,-1 sum=x(i) do 13 k=i+1,n sum=sum-a(k,i)*x(k) 13 continue x(i)=sum/p(i) 14 continue return END