program errorfldscoil ! vec pot, bx, by and bs for single coil. implicit double precision (a-h,o-z) c dimension vector(501,201),brfld(501,201),bsfld(501,201), + b1(501,201,360), b2(501,201,360), b3(501,201,360) dimension rad(201),xlong(501),bx(501,201),by(501,201), + bs(501,201) common/parameters/XL,R1 ! XL and R1 length of the field region c eta is position of center and rho radius external a1,a2,a3,a4,a5,a6,a7,a8 open(unit=29,file='pot_sol_coil.dat',status='unknown', + form='unformatted') open(unit=30,file='bx_sol_coil.dat',status='unknown', + form='unformatted') open(unit=31,file='by_sol_coil.dat',status='unknown', + form='unformatted') open(unit=32,file='bs_sol_coil.dat',status='unknown', + 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 coil in [m]' read(5,*)R1 c R1=0.7 XL=2*R1 ! total length is 2*XL; XL is 2 times the radius n_radial_points=200 n_long_points=500 ds=2*XL/real(n_long_points) dr=0.8*R1/real(n_radial_points) ! 80 % inside solenoid delta=0 ! radial position of coil ( i the paper is R) eta=0 ! long position of the coil (in the paper is S) c potential as a power series expansion in r c b00=mu_o*I/2 b00=0.5*perm_mu0*curr c order of the expansion c type*,'up to which order in r? (1,3,5,7...(le.7))' c read*,iorder write(6,*)'ENTER:tilt angle(psi/rad) and dir.of rotation(degree)' read(5,*)psi,roteta roteta=roteta*piover2/90.d0 theta=0. c write(6,*)'do you want to ENTER number of angles? Y(1);n(0)' c read(5,*)iloop c if(iloop.eq.1)then c write(6,*)'ENTER number of angles' c read(5,*)iangle write(6,*)'ENTER angle number for output (< 200)' c read(5,*)iout c else c iangle=360 c iout=1 c endif iangle=200 dtheta=1.d0/real(iangle) ! arbitrarily chosen do m=1,iangle+1 write(6,*)'angle value=',m, theta do k=1,n_radial_points+1 do i=1,n_long_points+1 vector(i,k)=0 brfld(i,k)=0 bsfld(i,k)=0 end do end do c loop over positions r=0. do k=1,n_radial_points+1 rad(k)=r s=-XL do i=1,n_long_points+1 xlong(i)=s c Now the rotated quantities x=r*cos(theta) y=r*sin(theta) aux1=x*cos(roteta)+y*sin(roteta) aux2=x*sin(roteta)-y*cos(roteta) pris=s*cos(psi)+aux2*sin(psi) c prir=sqrt(x**2+(y*cos(psi))**2+(s*sin(psi))**2- c + 2d0*s*y*cos(psi)*sin(psi)) prir=sqrt((r*cos(psi))**2-(aux1*(1d0-cos(psi)))**2+ + (s*sin(psi))**2+2.*cos(psi)*(1d0-cos(psi))*aux1**2- + 2.*cos(psi)*sin(psi)*s*aux2) c vec1=0.5*a1(s,eta) vec3=1d0/16.*a3(s,eta) vec5=1d0/384.*a5(s,eta) vec7=1d0/18432.*a7(s,eta) c br1=0.5*a2(s,eta) br2=1d0/16.*a4(s,eta) br3=1d0/384.*a6(s,eta) br4=1d0/18432.*a8(s,eta) c pribr1=0.5*a2(pris,eta) pribr2=1d0/16.*a4(pris,eta) pribr3=1d0/384.*a6(pris,eta) pribr4=1d0/18432.*a8(pris,eta) c bs1=a1(s,eta) bs2=1d0/8.*a3(s,eta) bs3=1d0/192.*a5(s,eta) bs4=1d0/9216.*a7(s,eta) c pribs1=a1(pris,eta) pribs2=1d0/8.*a3(pris,eta) pribs3=1d0/192.*a5(pris,eta) pribs4=1d0/9216.*a7(pris,eta) c c vector(i,k)=b00*r*(vec1+r**2*(vec3+r**2*(vec5+r**2*vec7))) c brfld(i,k)=-b00*r*(br1+r**2*(br2+r**2*(br3+r**2*br4))) c bsfld(i,k)=b00*(bs1+r**2*(bs2+r**2*(bs3+r**2*bs4))) c vector(i,k)=r*(vec1+r**2*(vec3+r**2*(vec5+r**2*vec7))) brfld(i,k)=-r*(br1+r**2*(br2+r**2*(br3+r**2*br4))) bsfld(i,k)=(bs1+r**2*(bs2+r**2*(bs3+r**2*bs4))) c c privector(i,k)=prir*(privec1+prir**2*(privec3+prir**2* c + (privec5+prir**2*privec7))) pribrfld=-prir*(pribr1+prir**2*(pribr2+prir**2* + (pribr3+prir**2*pribr4))) pribsfld=(pribs1+prir**2*(pribs2+prir**2* + (pribs3+prir**2*pribs4))) c rotated fields - the unrotated fields === ERROR FIELDS b1(i,k,m)=(pribrfld-brfld(i,k))*cos(roteta)+ + pribrfld*sin(psi)*sin(roteta)-pribsfld*sin(psi) b2(i,k,m)=(pribrfld*cos(psi)-brfld(i,k))*cos(roteta)- + pribrfld*sin(roteta)+sin(psi)*pribsfld b3(i,k,m)=-sin(psi)*pribrfld+ pribsfld*cos(psi)-bsfld(i,k) c b1(i,k,m)=pribrfld*sin(theta) c b2(i,k,m)=pribrfld*cos(theta)+psi*pribsfld c b3(i,k,m)=-psi*pribrfld*cos(theta)+ pribsfld c if(i.eq.300.and.m.eq.iout)then write(20,*)rad(k),b1(i,k,m) write(21,*)rad(k),b2(i,k,m) write(22,*)rad(k),b3(i,k,m) endif s=s+ds end do ! end s do loop if(k.eq.1.and.m.eq.iout)then do i=1,n_long_points+1 write(10,*)xlong(i),b1(i,k,m) write(11,*)xlong(i),b2(i,k,m) write(12,*)xlong(i),b3(i,k,m) end do endif r=r+dr end do ! end r do loop theta=theta+dtheta*4.*piover2 end do ! end theta do loop write(29)n_radial_points+1,n_long_points+1 write(29)vector write(30)n_radial_points+1,n_long_points+1 do m=1,iangle if(m.eq.iout)then do k=1,n_radial_points+1 do i=1,n_long_points+1 bx(i,k)=b1(i,k,m) end do end do write(30)bx endif end do write(31)n_radial_points+1,n_long_points+1 do m=1,iangle if(m.eq.iout)then do k=1,n_radial_points+1 do i=1,n_long_points+1 by(i,k)=b2(i,k,m) end do end do write(31)by endif end do write(32)n_radial_points+1,n_long_points+1 do m=1,iangle if(m.eq.iout)then do k=1,n_radial_points+1 do i=1,n_long_points+1 bs(i,k)=b3(i,k,m) end do end do write(32)bs endif end do write(6,*)'pot_sol_coil.dat contains:POT; NORMALIZE', + n_radial_points+1,n_long_points+1 write(6,*)'bx_sol_coil.dat contains:Bx; NORMALIZE', + n_radial_points+1,n_long_points+1 write(6,*)'by_sol_coil.dat contains:By; NORMALIZE', + n_radial_points+1,n_long_points+1 write(6,*)'bs_sol_coil.dat contains:Bs; NORMALIZE', + n_radial_points+1,n_long_points+1 write(6,*)'fort.10 contains: bx(1,s) vs s' write(6,*)'fort.11 contains: by(1,s) vs s' write(6,*)'fort.12 contains: bs(1,s) vs s' write(6,*)'fort.20 contains: bx(r,250) vs r' write(6,*)'fort.21 contains: by(r,250) vs r' write(6,*)'fort.22 contains: bs(r,250) vs r' c 50 format(1x,6(g12.5,1x)) stop end c FUNCTIONS real*8 function a1(ss,eta) ! 1st order implicit double precision (a-h,o-z) common/parameters/XL,R1 c eta is position of center and R1 radius aux0=sqrt(R1**2+(ss-eta)**2) aux1=1/aux0 aux2=-(ss-eta)**2/aux0**3 a1=aux1+aux2 end c real*8 function a2(ss,eta) ! 2nd order implicit double precision (a-h,o-z) common/parameters/XL,R1 c eta is position of center and R1 radius aux0=sqrt(R1**2+(ss-eta)**2) aux1=-3d0*(ss-eta)/aux0**3 aux2=3d0*(ss-eta)**3/aux0**5 a2=aux1+aux2 end c real*8 function a3(ss,eta) ! 3rd order implicit double precision (a-h,o-z) c eta is position of center and R1 radius common/parameters/XL,R1 aux0=sqrt(R1**2+(ss-eta)**2) aux1=-3d0/aux0**3 aux2=18d0*(ss-eta)**2/aux0**5 aux3=-15d0*(ss-eta)**4/aux0**7 a3=aux1+aux2+aux3 end c real*8 function a4(ss,eta) ! 4th order implicit double precision (a-h,o-z) common/parameters/XL,R1 c eta is position of center and R1 radius aux0=sqrt(R1**2+(ss-eta)**2) aux1=45d0*(ss-eta)/aux0**5 aux2=-150d0*(ss-eta)**3/aux0**7 aux3=105d0*(ss-eta)**5/aux0**9 a4=aux1+aux2+aux3 end c real*8 function a5(ss,eta) ! 5th order implicit double precision (a-h,o-z) c eta is position of center and R1 radius common/parameters/XL,R1 aux0=sqrt(R1**2+(ss-eta)**2) aux1=-945d0*(ss-eta)**6/aux0**11 aux2=1575d0*(ss-eta)**4/aux0**9 aux3=-675d0*(ss-eta)**2/aux0**7 aux4=45d0/aux0**5 a5=aux1+aux2+aux3+aux4 end c real*8 function a6(ss,eta) ! 6th order implicit double precision (a-h,o-z) common/parameters/XL,R1 c eta is position of center and R1 radius aux0=sqrt(R1**2+(ss-eta)**2) aux1=-1575d0*(ss-eta)/aux0**7 aux2=11025d0*(ss-eta)**3/aux0**9 aux3=-19845d0*(ss-eta)**5/aux0**11 aux4=10395d0*(ss-eta)**7/aux0**13 a6=aux1+aux2+aux3+aux4 end c real*8 function a7(ss,eta) ! 7th order implicit double precision (a-h,o-z) c eta is position of center and R1 radius common/parameters/XL,R1 aux0=sqrt(R1**2+(ss-eta)**2) aux1=-135135*(ss-eta)**8/aux0**15 aux2=291060*(ss-eta)**6/aux0**13 aux3=-198450*(ss-eta)**4/aux0**11 aux4=44100*(ss-eta)**2/aux0**9 aux5=-1575/aux0**7 a7=aux1+aux2+aux3+aux4+aux5 end c real*8 function a8(ss,eta) ! 8th order implicit double precision (a-h,o-z) common/parameters/XL,R1 c eta is position of center and R1 radius aux0=sqrt(R1**2+(ss-eta)**2) aux1=2027025d0*(ss-eta)**9/aux0**17 aux2=-4864860d0*(ss-eta)**7/aux0**15 aux3=3929310d0*(ss-eta)**5/aux0**13 aux4=-1190700d0*(ss-eta)**3/aux0**11 aux5=99225d0*(ss-eta)/aux0**9 a8=aux1+aux2+aux3+aux4+aux5 end c