program errorfldscoil ! 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 vector(201,501),rad(201),xlong(501) dimension b1(201,501),b2(201,501),b3(201,501) dimension bx(201,501),by(201,501),bs(201,501) dimension c1(201,501),c2(201,501),c3(201,501) dimension d1(201,501),d2(201,501),d3(201,501) dimension q1(201,501),q2(201,501),q3(201,501) dimension cx(201,501),cy(201,501),cs(201,501) dimension fx(201,501),fy(201,501),fs(201,501) dimension qx(201,501),qy(201,501),qs(201,501) 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') open(unit=60,file='cx_sol_coil.dat',status='unknown', + form='unformatted') open(unit=61,file='cy_sol_coil.dat',status='unknown', + form='unformatted') open(unit=62,file='cs_sol_coil.dat',status='unknown', + form='unformatted') open(unit=90,file='fx_sol_coil.dat',status='unknown', + form='unformatted') open(unit=91,file='fy_sol_coil.dat',status='unknown', + form='unformatted') open(unit=92,file='fs_sol_coil.dat',status='unknown', + form='unformatted') open(unit=10,file='qx_sol_coil.dat',status='unknown', + form='unformatted') open(unit=11,file='qy_sol_coil.dat',status='unknown', + form='unformatted') open(unit=12,file='qs_sol_coil.dat',status='unknown', + form='unformatted') c eta is position of center and rho radius away from axis 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 write(6,*)'Enter: length of region in [m]' read(5,*)slong c R1=0.7 XL=slong n_radial_points=200 n_long_points=500 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 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) & dir.of rotation(degree)' read(5,*)psi,roteta roteta=roteta*piover2/90.d0 write(6,*)'ENTER angle number for output (< 360)' read(5,*)iout theta=0. iangle=360 dtheta=4.*piover2/real(iangle) ! arbitrarily chosen c NOW translation of the coil write(6,*)'ENTER:trans. parameter[m] (angle: dir of rot.above)' read(5,*)trasla c define auxiliary functions trig1=sin(psi) trig2=cos(psi) trig3=sin(roteta) trig4=cos(roteta) c do m=1,iangle do k=1,n_radial_points+1 do i=1,n_long_points+1 b1(k,i)=0 b2(k,i)=0 b3(k,i)=0 bx(k,i)=0 by(k,i)=0 bs(k,i)=0 end do end do end do do m=1,iangle ! loop over field position angle if(mod(m,90).eq.0)write(6,*)'angle value=',m, theta c loop over positions r=0. do k=1,n_radial_points+1 ! radial position rad(k)=r s=-XL do i=1,n_long_points+1 ! longitudinal position xlong(i)=s c Now the rotated quantities x=r*cos(theta) y=r*sin(theta) aux1=x*trig4+y*trig3 aux2=x*trig3-y*trig4 pris=s*trig2+aux2*trig1 prix=x*trig2+trig4*aux1*(1-trig2)-s*trig1* + trig3 priy=y*trig2+trig3*aux1*(1-trig2)+s*trig2* + trig3 if(prix.eq.0.and.priy.eq.0)then prir=0 else prir=sqrt(prix**2+priy**2) endif c Now the translated quantities trax=x-trasla*trig4 tray=y-trasla*trig3 trar=sqrt(r**2+trasla**2-2.*trasla*(x*trig4+y*trig3)) tras=s c vec1=0.5*a1(s,eta) ! eta is S and delta is R vec3=-1.d0/16.*a3(s,eta) vec5=1.d0/384.*a5(s,eta) vec7=-1.d0/18432.*a7(s,eta) c NORMAL br1=0.5*a2(s,eta) br2=-1.d0/16.*a4(s,eta) br3=1.d0/384.*a6(s,eta) br4=-1.d0/18432.*a8(s,eta) c ROTATED pribr1=0.5*a2(pris,eta) pribr2=-1.d0/16.*a4(pris,eta) pribr3=1.d0/384.*a6(pris,eta) pribr4=-1.d0/18432.*a8(pris,eta) c TRANSLATED c no needed because s do not change c NORMAL bs1= a1(s,eta) bs2=-1.d0/4.*a3(s,eta) bs3=1d0/64.*a5(s,eta) bs4=-1d0/2304.*a7(s,eta) c ROTATED pribs1= a1(pris,eta) pribs2=-1.d0/4.*a3(pris,eta) pribs3=1.d0/64.*a5(pris,eta) pribs4=-1.d0/2304.*a7(pris,eta) c TRANSLATED c no needed because s do not change c c vector(k,i)=b00*r*(vec1+r**2*(vec3+r**2*(vec5+r**2*vec7))) c brfld(k,i)=-b00*r*(br1+r**2*(br2+r**2*(br3+r**2*br4))) c bsfld(k,i)=b00*(bs1+r**2*(bs2+r**2*(bs3+r**2*bs4))) c NORMAL vector(k,i)=b00*r*(vec1+r**2*(vec3+r**2*(vec5+r**2*vec7))) brfld=-b00*r*(br1+r**2*(br2+r**2*(br3+r**2*br4))) bsfld=b00*(bs1+r**2*(bs2+r**2*(bs3+r**2*bs4))) if(r.eq.0)then bxfld=0 byfld=0 else bxfld=brfld*x/r byfld=brfld*y/r endif cROTATION pribrfld=-b00*prir*(pribr1+prir**2*(pribr2+prir**2* + (pribr3+prir**2*pribr4))) pribsfld=b00*(pribs1+prir**2*(pribs2+prir**2* + (pribs3+prir**2*pribs4))) if(prir.eq.0)then pribxfld=0 pribyfld=0 else pribxfld=pribrfld*prix/prir pribyfld=pribrfld*priy/prir endif c Translation trabxfld= + -b00*trax*(br1+trar**2*(br2+trar**2*(br3+trar**2*br4))) trabyfld= + -b00*tray*(br1+trar**2*(br2+trar**2*(br3+trar**2*br4))) trabsfld= + b00*(bs1+trar**2*(bs2+trar**2*(bs3+trar**2*bs4))) c Rotation rotxfld=pribxfld*(trig2+trig4**2*(1.-trig2))+ + pribyfld*trig3*trig4*(1.-trig2)+ + pribsfld*trig3*trig1 rotyfld=pribxfld*trig3*trig4*(1.-trig2)+ + pribyfld*(trig2+trig3**2*(1.-trig2))+ + pribsfld*trig4*trig1 rotsfld=trig1*(pribxfld*trig3- + pribyfld*trig4)+ + pribsfld*trig2 c rotated fields - the unrotated fields === ERROR FIELDS b1(k,i)= rotxfld- pribxfld b2(k,i)= rotyfld- pribyfld b3(k,i)= rotsfld- pribsfld c ROTATED FIELDS c1(k,i)= rotxfld c2(k,i)= rotyfld c3(k,i)= rotsfld c UNROTATED FIELDS d1(k,i)= pribxfld d2(k,i)= pribyfld d3(k,i)= pribsfld c ERROR field due to translation q1(k,i)=bxfld-trabxfld q2(k,i)=byfld-trabyfld q3(k,i)=bsfld-trabsfld c if(i.eq.imiddle.and.m.eq.iout)then write(20,*)rad(k),b1(k,i),q1(k,i) write(21,*)rad(k),b2(k,i),q2(k,i) write(22,*)rad(k),b3(k,i),q3(k,i) endif 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(23,*)xlong(i),b1(k,i),q1(k,i) write(24,*)xlong(i),b2(k,i),q2(k,i) write(25,*)xlong(i),b3(k,i),q3(k,i) end do endif r=r+dr end do ! end r do loop theta=theta+dtheta if(m.eq.iout)then do k=1,n_radial_points+1 do i=1,n_long_points+1 bx(k,i)=b1(k,i) by(k,i)=b2(k,i) bs(k,i)=b3(k,i) qx(k,i)=q1(k,i) qy(k,i)=q2(k,i) qs(k,i)=q3(k,i) cx(k,i)=c1(k,i) cy(k,i)=c2(k,i) cs(k,i)=c3(k,i) fx(k,i)=d1(k,i) fy(k,i)=d2(k,i) fs(k,i)=d3(k,i) end do end do endif end do ! end theta do loop c write(29)n_radial_points+1,n_long_points+1 write(29)vector ! it does not depend on angle theta write(30)n_radial_points+1,n_long_points+1 write(30)bx write(60)cx write(90)fx write(10)qx write(31)n_radial_points+1,n_long_points+1 write(31)by write(61)cy write(91)fy write(11)qy write(32)n_radial_points+1,n_long_points+1 write(32)bs write(62)cs write(92)fs write(12)qs write(6,*)'pot_sol_coil.dat contains:POT; NORMALIZED', + n_radial_points+1,n_long_points+1 write(6,*)'bx_sol_coil.dat contains:Bx; NORMALIZED ROT-ERROR', + n_radial_points+1,n_long_points+1 write(6,*)'by_sol_coil.dat contains:By; NORMALIZED ROT-ERROR', + n_radial_points+1,n_long_points+1 write(6,*)'bs_sol_coil.dat contains:Bs; NORMALIZED ROT-ERROR', + n_radial_points+1,n_long_points+1 write(6,*)'qx_sol_coil.dat contains:Bx; NORMALIZED TRAS-ERROR', + n_radial_points+1,n_long_points+1 write(6,*)'qy_sol_coil.dat contains:By; NORMALIZED TRAS-ERROR', + n_radial_points+1,n_long_points+1 write(6,*)'qs_sol_coil.dat contains:Bs; NORMALIZED TRAS-ERROR', + n_radial_points+1,n_long_points+1 write(6,*)'cx_sol_coil.dat contains:Bx; NORMALIZED ROT', + n_radial_points+1,n_long_points+1 write(6,*)'cy_sol_coil.dat contains:By; NORMALIZED ROT', + n_radial_points+1,n_long_points+1 write(6,*)'cs_sol_coil.dat contains:Bs; NORMALIZED ROT', + n_radial_points+1,n_long_points+1 write(6,*)'fx_sol_coil.dat contains:Bx; NORMALIZED UNROT', + n_radial_points+1,n_long_points+1 write(6,*)'fy_sol_coil.dat contains:By; NORMALIZED UNROT', + n_radial_points+1,n_long_points+1 write(6,*)'fs_sol_coil.dat contains:Bs; NORMALIZED UNROT', + n_radial_points+1,n_long_points+1 c 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