program helicalfld ! vec pot, bx, by and bs for single coil. implicit double precision (a-h,o-z) common/parameters/R1,alpha,f0 !XL and R1 length of the field region common/coord/r,theta,s c external ff1,ff2,ff3 dimension rad(181),xlong(181) dimension b1(181,181,181),b2(181,181,181),b3(181,181,181) dimension bxfld(181,181),byfld(181,181),bsfld(181,181) dimension brfld(181,181),bffld(181,181) c open(unit=27,file='br_helix.dat',status='unknown', + form='unformatted') open(unit=28,file='bf_helix.dat',status='unknown', + form='unformatted') open(unit=29,file='bs_helix.dat',status='unknown', + form='unformatted') piover2=asin(1.0d0) open(unit=30,file='bx_helix.dat',status='unknown', + form='unformatted') open(unit=31,file='by_helix.dat',status='unknown', + form='unformatted') 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 helical coil in [m]' read(5,*)xlength write(6,*)'Enter: length of region in [m]' read(5,*)slong c alpha=4*piover2/xlambda alpha_radius=alpha*R1 XL=slong f0=0. ! initial phase in other words different helixs write(6,*)'Enter number of points: radial,long, angle' read(5,*)n_radial_points,n_long_points, n_angle c n_radial_pointsn_radial_points=180 c n_long_points=180 c n_angle=180 imiddle=n_long_points/2 kmiddle=n_radial_points/2 jmiddle=n_angle/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 dx=0.5*xlength c perm_mu0=8*piover2*1.d-7 ! permeability c b00=perm_mu0*curr/(8*piover2) b00=1.d-7*curr x00=b00*2*alpha_radius c write(6,*)'ENTER angle number for output (< 360)' c read(5,*)iout theta=0.d0 c loop field position angle c n_angle=1 do m=1,n_angle if(mod(m,jmiddle).eq.0)write(6,*)'angle value=',m, theta cost=cos(theta) sint=sin(theta) r=0.d0 c loop over positions do k=1,n_radial_points+1 ! radial position rad(k)=r c x=r*cost c y=r*sint s=-XL do i=1,n_long_points+1 ! longitudinal position xlong(i)=s space=-xlength ! length of the helical coil resul1=0 resul2=0 resul3=0 do while(space.le.xlength) spaceend=space+dx call qromb(ff1,space,spaceend,result1) call qromb(ff2,space,spaceend,result2) call qromb(ff3,space,spaceend,result3) c write(6,*)k,i,result3 resul1=resul1+result1 resul2=resul2+result2 resul3=resul3+result3 space=spaceend end do b1(k,m,i)=x00/alpha_radius*resul1 b2(k,m,i)=x00/alpha_radius*resul2 b3(k,m,i)=x00*resul3 c if(m.eq.jmiddle)then brfld(k,i)=b1(k,m,i) bffld(k,i)=b2(k,m,i) bsfld(k,i)=b3(k,m,i) bxfld(k,i)=-b1(k,m,i)*sint+b2(k,m,i)*cost byfld(k,i)=b1(k,m,i)*cost+b2(k,m,i)*sint if(i.eq.1)then write(20,*)rad(k),b1(k,m,i) write(21,*)rad(k),b2(k,m,i) write(22,*)rad(k),b3(k,m,i) else if(i.eq.imiddle)then write(20,*)rad(k),b1(k,m,i) write(21,*)rad(k),b2(k,m,i) write(22,*)rad(k),b3(k,m,i) else if(i.eq.n_long_points)then write(20,*)rad(k),b1(k,m,i) write(21,*)rad(k),b2(k,m,i) write(22,*)rad(k),b3(k,m,i) endif endif s=s+ds end do ! end s do loop if(m.eq.jmiddle)then if(k.eq.1)then do i=1,n_long_points+1 write(24,*)xlong(i),b1(k,m,i) write(25,*)xlong(i),b2(k,m,i) write(26,*)xlong(i),b3(k,m,i) end do else if(k.eq.kmiddle)then do i=1,n_long_points+1 write(24,*)xlong(i),b1(k,m,i) write(25,*)xlong(i),b2(k,m,i) write(26,*)xlong(i),b3(k,m,i) end do else if(k.eq.n_radial_points)then do i=1,n_long_points+1 write(24,*)xlong(i),b1(k,m,i) write(25,*)xlong(i),b2(k,m,i) write(26,*)xlong(i),b3(k,m,i) end do endif endif r=r+dr end do ! end r do loop theta=theta+dtheta end do ! end theta do loop c write(27)bxfld write(28)byfld write(29)bsfld write(30)bxfld write(31)byfld c write(6,*)'fort.20 contains: b1(r,iout,100) vs r' write(6,*)'fort.21 contains: b2(r,iout,100) vs r' write(6,*)'fort.22 contains: b3(r,iout,100) vs r' write(6,*)'fort.24 contains: b1(1,iout,s) vs s' write(6,*)'fort.25 contains: b2(100,iout,s) vs s' write(6,*)'fort.26 contains: b3(100,iout,s) vs s' write(6,*)'br_helix.dat contains:A(i,j)', + n_radial_points+1,n_long_points+1 write(6,*)'bf_helix.dat contains:A(i,j)', + n_radial_points+1,n_long_points+1 write(6,*)'bs_helix.dat contains:A(i,j)', + n_radial_points+1,n_long_points+1 write(6,*)'bx_helix.dat contains:A(i,j)', + n_radial_points+1,n_long_points+1 write(6,*)'by_helix.dat contains:A(i,j)', + n_radial_points+1,n_long_points+1 stop end c real*8 function ff1(t) implicit double precision (a-h,o-z) common/parameters/R1,alpha,f0!XL and R1 length field region common/coord/r,theta,s xyz=cos(theta-alpha*t+f0) aux2=r*r+R1*(R1-2.*r*xyz)+(s-t)**2 aux1=alpha*(s-t)*xyz-sin(theta-alpha*t+f0) ff1=R1*aux1/aux2**1.5 end c real*8 function ff2(t) implicit double precision (a-h,o-z) common/parameters/R1,alpha,f0 !XL and R1 length field region common/coord/r,theta,s xyz=cos(theta-alpha*t+f0) aux1=r-R1*(xyz+alpha*(s-t)*sin(theta-alpha*t+f0)) aux2=r*r+R1*(R1-2.*r*xyz)+(s-t)**2 ff2= aux1/aux2**1.5 end c function ff3(t) implicit double precision (a-h,o-z) common/parameters/R1,alpha,f0 ! XL and R1 length field region common/coord/r,theta,s xyz=cos(theta-alpha*t+f0) aux1=-r*xyz+R1 aux2=r*r+R1*(R1-2.*r*xyz)+(s-t)**2 ff3=aux1/aux2**1.5 end c SUBROUTINE polint(xa,ya,n,x,y,dy) 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=abs(x-xa(1)) do 11 i=1,n dift=abs(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.d0)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 C (C) Copr. 1986-92 Numerical Recipes Software *#'".,. c SUBROUTINE qromb(func,a,b,ss) INTEGER JMAX,JMAXP,K,KM DOUBLE PRECISION a,b,func,ss,EPS EXTERNAL func PARAMETER (EPS=1.d-6, JMAX=20, JMAXP=JMAX+1, K=5, KM=K-1) CU USES polint,trapzd INTEGER j DOUBLE PRECISION dss,h(JMAXP),s(JMAXP) h(1)=1.d0 do 11 j=1,JMAX call trapzd(func,a,b,s(j),j) if (j.ge.K) then call polint(h(j-KM),s(j-KM),K,0.d0,ss,dss) if (abs(dss).le.EPS*abs(ss)) return endif s(j+1)=s(j) h(j+1)=0.25d0*h(j) 11 continue pause 'too many steps in qromb' END C (C) Copr. 1986-92 Numerical Recipes Software *#'".,. c SUBROUTINE trapzd(func,a,b,s,n) INTEGER n DOUBLE PRECISION a,b,s,func EXTERNAL func INTEGER it,j DOUBLE PRECISION del,sum,tnm,x if (n.eq.1) then s=0.5d0*(b-a)*(func(a)+func(b)) else it=2**(n-2) tnm=it del=(b-a)/tnm x=a+0.5d0*del sum=0.d0 do 11 j=1,it sum=sum+func(x) x=x+del 11 continue s=0.5d0*(s+(b-a)*sum/tnm) endif return END C (C) Copr. 1986-92 Numerical Recipes Software *#'".,.