program helicalfld ! vec pot, bx, by and bs for single coil. implicit double precision (a-h,o-z) common/parameters/XL,R1 ! XL and R1 length of the field region c external a1,a2,a3,a4,a5,a6,a7,a8 dimension scalar(361,361),rad(361),xlong(361) dimension brfld(361,361),bphifld(361,361),bsfld(361,361) dimension b1(361,361,361),b2(361,361,361),b3(361,361,361) c open(unit=12,file='qs_sol_coil.dat',status='unknown', c + form='unformatted') piover2=asin(1.0d0) perm_mu0=8*piover2*1.d-7 ! permeability write(6,*)'enter: current in [A]' read(5,*)curr write(6,*)'Enter: radius of helical coil in [m]' read(5,*)R1 write(6,*)'Enter: period of helical coil in [m]' read(5,*)xlambda write(6,*)'Enter: length of region in [m]' read(5,*)slong c alpha=4*piover2/xlambda alpha_radius=alpha*R1 XL=slong n_radial_points=360 n_long_points=360 n_angle=360 imiddle=n_long_points/2 kmiddle=n_radial_points/2 ds=2*XL/real(n_long_points) dr=0.8*R1/real(n_radial_points) ! 80 % inside solenoid dtheta=4.*piover2/real(n_angle) ! arbitrarily chosen b00=2.*perm_mu0*curr x00=b00*2*alpha_radius c potential as a power series expansion in r c order of the expansion write(6,*)'up to which order in r? (1,3,5,7...(le.7))' read(5,*)korder c write(6,*)'ENTER:tilt angle(psi/rad) & dir.of rotation(degree)' c read(5,*)psi,roteta c roteta=roteta*piover2/90.d0 c NOW translation of the coil c write(6,*)'ENTER:trans. parameter[m] (angle: dir of rot.above)' c read(5,*)trasla write(6,*)'ENTER angle number for output (< 360)' read(5,*)iout theta=0. c loop field position angle do m=1,n_angle if(mod(m,45).eq.0)write(6,*)'angle value=',m, theta cost=cos(theta) sint=sin(theta) c loop over positions r=0. do k=1,n_radial_points+1 ! radial position rad(k)=r x=r*cost y=r*sint s=-XL do i=1,n_long_points+1 ! longitudinal position xlong(i)=s c do iorder=1,korder aux1=iorder*alpha_radius aux2=iorder*alpha*r aux3=iorder*(theta-alpha*s) sinsin=sin(aux3) coscos=cos(aux3) if(iorder.eq.1)then besskder=-0.5*(bessk0(aux1)+bessk(iorder+1,aux1)) bessider= 0.5*(bessi0(aux2)+bessi(iorder+1,aux2)) if(r.eq.0)then besseli=alpha*0.5 else besseli=bessi1(aux2)/r endif else if(iorder.eq.2)then besskder=-0.5*(bessk1(aux1)+bessk(iorder+1,aux1)) bessider= 0.5*(bessi1(aux2)+bessi(iorder+1,aux2)) if(r.eq.0)then besseli=0. else besseli=bessi(iorder,aux2)/r endif else besskder=-0.5*(bessk(iorder-1,aux1)+bessk(iorder+1,aux1)) bessider= 0.5*(bessi(iorder-1,aux2)+bessi(iorder+1,aux2)) if(r.eq.0)then besseli=0. else besseli=bessi(iorder,aux2)/r endif endif c sca=sca+besskder*besseli*sinsin br= br+besskder*bessider*sinsin*iorder bphi=bphi+besskder*besseli*iorder*coscos bs=bs+besskder*besseli*iorder*alpha*coscos end do c scalar(k,i)=b00*alpha*s+x00*sca brfld(k,i)=x00*br*alpha bphifld(k,i)=x00*bphi bsfld(k,i)=b00*alpha-x00*bs c if(i.eq.imiddle.and.m.eq.iout)then write(20,*)rad(k),scalar(k,i) write(21,*)rad(k),brfld(k,i) write(22,*)rad(k),bphifld(k,i) write(23,*)rad(k),bsfld(k,i) endif b1(k,m,i)=brfld(k,i)*cost-bphifld(k,i)*sint b2(k,m,i)=brfld(k,i)*sint+bphifld(k,i)*cost b3(k,m,i)=bsfld(k,i) s=s+ds end do ! end s do loop if(k.eq.kmiddle.and.m.eq.iout)then do i=1,n_long_points+1 write(24,*)xlong(i),scalar(k,i) write(25,*)xlong(i),brfld(k,i) write(26,*)xlong(i),bphifld(k,i) write(27,*)xlong(i),bsfld(k,i) end do endif r=r+dr end do ! end r do loop theta=theta+dtheta c if(m.eq.iout)then c do k=1,n_radial_points+1 c do i=1,n_long_points+1 c bx(k,i)=b1(k,i) c by(k,i)=b2(k,i) c bs(k,i)=b3(k,i) c qx(k,i)=q1(k,i) c qy(k,i)=q2(k,i) c qs(k,i)=q3(k,i) c cx(k,i)=c1(k,i) c cy(k,i)=c2(k,i) c cs(k,i)=c3(k,i) c fx(k,i)=d1(k,i) c fy(k,i)=d2(k,i) c fs(k,i)=d3(k,i) c end do c end do c endif end do ! end theta do loop c c write(6,*)'fort.20 contains: potential(r,iout,100) vs r' write(6,*)'fort.21 contains: b1(r,iout,100) vs r' write(6,*)'fort.22 contains: b2(r,iout,100) vs r' write(6,*)'fort.23 contains: b3(r,iout,100) vs r' write(6,*)'fort.24 contains: potential(1,iout,s) vs s' write(6,*)'fort.25 contains: b1(100,iout,s) vs s' write(6,*)'fort.26 contains: b2(100,iout,s) vs s' write(6,*)'fort.27 contains: b3(100,iout,s) vs s' stop end FUNCTION bessk(n,x) INTEGER n DOUBLE PRECISION bessk,x CU USES bessk0,bessk1 INTEGER j DOUBLE PRECISION bk,bkm,bkp,tox,bessk0,bessk1 if (n.lt.2) pause 'bad argument n in bessk' tox=2.0d0/x bkm=bessk0(x) bk=bessk1(x) do 11 j=1,n-1 bkp=bkm+j*tox*bk bkm=bk bk=bkp 11 continue bessk=bk return END C (C) Copr. 1986-92 Numerical Recipes Software *#'".,. FUNCTION bessi(n,x) INTEGER n,IACC DOUBLE PRECISION bessi,x,BIGNO,BIGNI PARAMETER (IACC=40,BIGNO=1.0d10,BIGNI=1.0d-10) CU USES bessi0 INTEGER j,m DOUBLE PRECISION bi,bim,bip,tox,bessi0 if (n.lt.2) pause 'bad argument n in bessi' if (x.eq.0.d0) then bessi=0.d0 else tox=2.0d0/abs(x) bip=0.0d0 bi=1.0d0 bessi=0.d0 m=2*((n+int(sqrt( dble(IACC*n))))) do 11 j=m,1,-1 bim=bip+ dble(j)*tox*bi bip=bi bi=bim if (abs(bi).gt.BIGNO) then bessi=bessi*BIGNI bi=bi*BIGNI bip=bip*BIGNI endif if (j.eq.n) bessi=bip 11 continue bessi=bessi*bessi0(x)/bi if (x.lt.0.d0.and.mod(n,2).eq.1) bessi=-bessi endif return END C (C) Copr. 1986-92 Numerical Recipes Software *#'".,. FUNCTION bessi0(x) DOUBLE PRECISION bessi0,x DOUBLE PRECISION ax DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y * SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,3.5156229d0,3.0899424d0, *1.2067492d0,0.2659732d0,0.360768d-1,0.45813d-2/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,0.1328592d-1, *0.225319d-2,-0.157565d-2,0.916281d-2,-0.2057706d-1,0.2635537d-1, *-0.1647633d-1,0.392377d-2/ if (abs(x).lt.3.75d0) then y=(x/3.75d0)**2 bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))) else ax=abs(x) y=3.75d0/ax bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y* *(q7+y*(q8+y*q9)))))))) endif return END C (C) Copr. 1986-92 Numerical Recipes Software *#'".,. FUNCTION bessk0(x) DOUBLE PRECISION bessk0,x CU USES bessi0 DOUBLE PRECISION bessi0 DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7 DATA p1,p2,p3,p4,p5,p6,p7/-0.57721566d0,0.42278420d0,0.23069756d0, * *0.3488590d-1,0.262698d-2,0.10750d-3,0.74d-5/ DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,-0.7832358d-1,0.2189568d-1, * *-0.1062446d-1,0.587872d-2,-0.251540d-2,0.53208d-3/ if (x.le.2.0d0) then y=x*x/4.0d0 bessk0=(-log(x/2.0d0)*bessi0(x))+(p1+y*(p2+y*(p3+y*(p4+y*(p5+y* *(p6+y*p7)))))) else y=(2.0d0/x) bessk0=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y* *q7)))))) endif return END C (C) Copr. 1986-92 Numerical Recipes Software *#'".,. FUNCTION bessk1(x) DOUBLE PRECISION bessk1,x CU USES bessi1 DOUBLE PRECISION bessi1 DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7 DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,0.15443144d0,-0.67278579d0, *-0.18156897d0,-0.1919402d-1,-0.110404d-2,-0.4686d-4/ DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,0.23498619d0,-0.3655620d-1, * *0.1504268d-1,-0.780353d-2,0.325614d-2,-0.68245d-3/ if (x.le.2.0d0) then y=x*x/4.0d0 bessk1=(log(x/2.0d0)*bessi1(x))+(1.0d0/x)*(p1+y*(p2+y*(p3+y*(p4 *+y* *(p5+y*(p6+y*p7)))))) else y=2.0d0/x bessk1=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y* *q7)))))) endif return END C (C) Copr. 1986-92 Numerical Recipes Software *#'".,. FUNCTION bessi1(x) DOUBLE PRECISION bessi1,x DOUBLE PRECISION ax DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y * SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0, *0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,-0.3988024d-1, *-0.362018d-2,0.163801d-2,-0.1031555d-1,0.2282967d-1,-0.2895312d-1, * *0.1787654d-1,-0.420059d-2/ if (abs(x).lt.3.75d0) then y=(x/3.75d0)**2 bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else ax=abs(x) y=3.75d0/ax bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y* *(q7+y*(q8+y*q9)))))))) if(x.lt.0.d0)bessi1=-bessi1 endif return END C (C) Copr. 1986-92 Numerical Recipes Software *#'".,.