program vector_potential_sheet ! vector potential for thick c cylindrical current sheet. implicit double precision (a-h,o-z) c dimension vector(201,501), fldr(201,501), flds(201,501) dimension rad(201),xlong(501) common/parameters/XL,R1,R2 c R1 internal radius; R2 external radius external a1,a3,a5,a7 open(unit=29,file='pot_sheet.dat',status='unknown', + form='unformatted') open(unit=28,file='br_sheet.dat',status='unknown', + form='unformatted') open(unit=27,file='bs_sheet.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: radius1 of coil in [m]' read(5,*)R1 write(6,*)'Enter: radius2 of coil in [m]' read(5,*)R2 write(6,*)'Enter:length of magnet in [m](correspond to half' c write(6,*)'Enter: length of magnet in [m]' read(5,*)XL write(6,*)'Enter: length of region in [m]' read(5,*)slong write(6,*)'Enter: radius of region in [m]' read(5,*)rlong c n_radial_points=200 n_long_points=500 ds=2*slong/n_long_points dr=rlong/n_radial_points write(6,*)'dr= ',dr,' ds= ',ds imiddle=n_long_points/2 kmiddle=n_radial_points/2 eta=0 ! centered at origen c potential as a power series expansion in r c c X00=perm_mu0*curr/2./(R2-R1)*Log(R2+sqrt(R2**2+XL**2)/ c + R1+sqrt(R1**2+XL**2)) c b00=log((R2+sqrt(R2**2+XL**2))/(R1+sqrt(R1**2+XL**2)))/(2*XL) c b0=perm_mu0*curr/2. c b000=perm_mu0*curr/4./XL/(R2-R1) b00=perm_mu0*curr/4./XL/(R2-R1) !sign to have better graphics X00=2.*b00*XL*(Log(R2+sqrt(R2**2+XL**2))- + Log(R1+sqrt(R1**2+XL**2))) X01=2.*b00*XL*Log((R2+sqrt(R2**2+XL**2))/ + (R1+sqrt(R1**2+XL**2))) write(6,*)X00,'X00',X01 c both expression above are for cylindrical current sheet c b00=mu_o*I/2/R ! for coil c b00=R1 ! for coil c order of the expansion c type*,'up to which order in r? (1,3,5,7...(le.7))' c read*,iorder c collection of constants upot3=1.-(R1/R2)**3 upot4=1.-(R1/R2)**4 upot5=1.-(R1/R2)**5 upot6=1.-(R1/R2)**6 upot7=1.-(R1/R2)**7 upot8=1.-(R1/R2)**8 upot9=1.-(R1/R2)**9 upot10=1.-(R1/R2)**10 upot11=1.-(R1/R2)**11 do k=1,n_radial_points+1 do i=1,n_long_points+1 vector(k,i)=0 end do end do c loop over positions r=0 c if(r.eq.0)then vec1=(a1(0,R2,XL)-a1(0,R2,-XL)- + (a1(0,R1,XL)-a1(0,R1,-XL))) write(6,*)vec1,x00/b00, 'Normalized to Bs(r=0,s=0)=1 ' c b00=-b00/X00 b00=-b00 c endif do k=1,n_radial_points+1 rad(k)=r s=-slong c test for radial position if(r.le.R1)then do i=1,n_long_points+1 xlong(i)=s vec1=(a1(s,R2,XL)-a1(s,R2,-XL)- + (a1(s,R1,XL)-a1(s,R1,-XL))) c vec2=(a2(s,R2,XL)-a2(s,R2,-XL)- + (a2(s,R1,XL)-a2(s,R1,-XL))) c vec3=(a3(s,R2,XL)-a3(s,R2,-XL)- + (a3(s,R1,XL)-a3(s,R1,-XL)) ) c vec4=(a4(s,R2,XL)-a4(s,R2,-XL)- + (a4(s,R1,XL)-a4(s,R1,-XL)) ) c vec5=(a5(s,R2,XL)-a5(s,R2,-XL)- + (a5(s,R1,XL)-a5(s,R1,-XL))) c vec6=(a6(s,R2,XL)-a6(s,R2,-XL)- + (a6(s,R1,XL)-a6(s,R1,-XL))) c vec7=(a7(s,R2,XL)-a7(s,R2,-XL)- + (a7(s,R1,XL)-a7(s,R1,-XL))) vec8=(a8(s,R2,XL)-a8(s,R2,-XL)- + (a8(s,R1,XL)-a8(s,R1,-XL))) vec9=(a9(s,R2,XL)-a9(s,R2,-XL)- + (a9(s,R1,XL)-a9(s,R1,-XL))) c POTENTIAL vector(k,i)=-b00*r*(0.5*vec1+r**2*(-vec3/16.+ + r**2*(vec5/384.+r**2*(-vec7/18432.+r**2*vec9/5406720.)))) c RADIAL COMPONENT fldr(k,i)=b00*r*(0.5*vec2+r**2*(-vec4/16.+ + r**2*(vec6/384.-r**2*vec8/18432.))) c LONG COMPONENT flds(k,i)=-b00*(vec1+r**2*(-vec3/4.+r**2*(vec5/64.+ + r**2*(-vec7/2304.+r*2*vec9/540672.)))) if(i.eq.imiddle)then write(40,*)rad(k),vector(k,i) write(41,*)rad(k),fldr(k,i) write(42,*)rad(k),flds(k,i) endif s=s+ds end do if(k.eq.1)then do i=1,n_long_points+1 write(30,*)xlong(i),vector(k,i) write(31,*)xlong(i),fldr(k,i) write(32,*)xlong(i),flds(k,i) end do else if(k.eq.kmiddle)then do i=1,n_long_points+1 write(30,*)xlong(i),vector(k,i) write(31,*)xlong(i),fldr(k,i) write(32,*)xlong(i),flds(k,i) end do else if(k.eq.n_rad_points+1)then do i=1,n_long_points+1 write(30,*)xlong(i),vector(k,i) write(31,*)xlong(i),fldr(k,i) write(32,*)xlong(i),flds(k,i) end do endif else if (r.lt.R2) then ! if smaller than external radius R2 c s=-slong do i=1,n_long_points+1 xlong(i)=s tec1=q1(r,s,XL)-q1(r,s,-XL) uec1=z1(r,s,XL)-z1(r,s,-XL) pot3=1.-(R1/r)**3 vec1=(a1(s,R2,XL)-a1(s,R2,-XL)- + (a1(s,r,XL)-a1(s,r,-XL))) c tec2=q2(r,s,XL)-q2(r,s,-XL) uec2=z2(r,s,XL)-z2(r,s,-XL) pot4=1.-(R1/r)**4 vec2=(a2(s,R2,XL)-a2(s,R2,-XL)- + (a2(s,r,XL)-a2(s,r,-XL))) c tec3=q3(r,s,XL)-q3(r,s,-XL) uec3=z3(r,s,XL)-z3(r,s,-XL) pot5=1.-(R1/r)**5 vec3=(a3(s,R2,XL)-a3(s,R2,-XL)- + (a3(s,r,XL)-a3(s,r,-XL))) c tec4=q4(r,s,XL)-q4(r,s,-XL) uec4=z4(r,s,XL)-z4(r,s,-XL) pot6=1.-(R1/r)**6 vec4=(a4(s,R2,XL)-a4(s,R2,-XL)- + (a4(s,r,XL)-a4(s,r,-XL))) c tec5=q5(r,s,XL)-q5(r,s,-XL) uec5=z5(r,s,XL)-z5(r,s,-XL) pot7=1.-(R1/r)**7 vec5=(a5(s,R2,XL)-a5(s,R2,-XL)- + (a5(s,r,XL)-a5(s,r,-XL))) c tec6=q6(r,s,XL)-q6(r,s,-XL) uec6=z6(r,s,XL)-z6(r,s,-XL) pot8=1.-(R1/r)**8 vec6=(a6(s,R2,XL)-a6(s,R2,-XL)- + (a6(s,r,XL)-a6(s,r,-XL))) c tec7=q7(r,s,XL)-q7(r,s,-XL) uec7=z7(r,s,XL)-z7(r,s,-XL) pot9=1.-(R1/r)**9 vec7=(a7(s,R2,XL)-a7(s,R2,-XL)- + (a7(s,r,XL)-a7(s,r,-XL))) c tec8=q8(r,s,XL)-q8(r,s,-XL) uec8=z8(r,s,XL)-z8(r,s,-XL) pot10=1.-(R1/r)**10 vec8=(a8(s,R2,XL)-a8(s,R2,-XL)- + (a8(s,r,XL)-a8(s,r,-XL))) c tec9=q9(r,s,XL)-q9(r,s,-XL) uec9=z9(r,s,XL)-z9(r,s,-XL) pot11=1.-(R1/r)**11 vec9=(a9(s,R2,XL)-a9(s,R2,-XL)- + (a9(s,r,XL)-a9(s,r,-XL))) c VECTOR POTENTIAL c vector(k,i)=b00*(R2-R1)/(r-R1)* vector(k,i)=b00* + r**2*(tec1*pot3/6.+r**2*(-tec3*pot5/80.+ + r**2*(tec5*pot7/2688.+r**2*(-tec7*pot9/165888. + +r**2*tec9*pot11/8110080.)))) ! r larger than R1 c + +b00*(R2-R1)/(R2-r)*r*(0.5*vec1+r**2*(-vec3/16.+ + -b00*r*(0.5*vec1+r**2*(-vec3/16.+ + r**2*(vec5/384.+r**2*(-vec7/18432.+r**2*vec9/5406720.)))) c there is factor of 1/r in the functions tec* which is explicitly taken out c RADIAL MAGNETIC FIELD c fldr(k,i)=-b00*(R2-R1)/(1.-R1/r)* fldr(k,i)=b00* ! + r**2*(tec2*pot3/6.+r**2*(-tec4*pot5/80.+ + r**2*(tec6*pot7/2688.+r**2*(-tec8*pot9/165888.)))) c + +b00*(R2-R1)/(R2-r)*r*(0.5*vec2+r**2*(-vec4/16.+ + +b00*r*(0.5*vec2+r**2*(-vec4/16.+ + r**2*(vec6/384.-r**2*vec8/18432.))) c LONG MAGNETIC FIELD c flds(k,i)=b00*(R2-R1)/(1.-R1/r)* flds(k,i)=-b00* + r**3*(pot3*(uec1)/6.+r**2*(-pot5*(uec3)/80.+ + r**2*(pot7*(uec5)/2688.+r**2*(-pot9*(uec7)/165888. + +r**2*pot11*(uec9)/8110080.)))) c + -b00*(R2-R1)/(R2-r)*(vec1+r**2*(-vec3/4.+ + -b00*(vec1+r**2*(-vec3/4.+ + r**2*(vec5/64.+r**2*(-vec7/2304.+r**2*vec9/540672.)))) c if(i.eq.imiddle)then write(40,*)rad(k),vector(k,i) write(41,*)rad(k),fldr(k,i) write(42,*)rad(k),flds(k,i) endif s=s+ds end do if(k.eq.1)then do i=1,n_long_points+1 write(30,*)xlong(i),vector(k,i) write(31,*)xlong(i),fldr(k,i) write(32,*)xlong(i),flds(k,i) end do else if(k.eq.kmiddle)then do i=1,n_long_points+1 write(30,*)xlong(i),vector(k,i) write(31,*)xlong(i),fldr(k,i) write(32,*)xlong(i),flds(k,i) end do else if(k.eq.n_radial_points+1)then do i=1,n_long_points+1 write(30,*)xlong(i),vector(k,i) write(31,*)xlong(i),fldr(k,i) write(32,*)xlong(i),flds(k,i) end do endif else ! if r >= R2 do i=1,n_long_points+1 xlong(i)=s tec1=q1(r,s,XL)-q1(r,s,-XL) uec1=z1(r,s,XL)-z1(r,s,-XL) c upot3=1.-(R1/R2)**3 c tec2=q2(r,s,XL)-q2(r,s,-XL) uec2=z2(r,s,XL)-z2(r,s,-XL) c upot4=1.-(R1/R2)**4 c tec3=q3(r,s,XL)-q3(r,s,-XL) uec3=z3(r,s,XL)-z3(r,s,-XL) c upot5=1.-(R1/R2)**5 c tec4=q4(r,s,XL)-q4(r,s,-XL) uec4=z4(r,s,XL)-z4(r,s,-XL) c upot6=1.-(R1/R2)**6 c tec5=q5(r,s,XL)-q5(r,s,-XL) uec5=z5(r,s,XL)-z5(r,s,-XL) c upot7=1.-(R1/R2)**7 c tec6=q6(r,s,XL)-q6(r,s,-XL) uec6=z6(r,s,XL)-z6(r,s,-XL) c upot8=1.-(R1/R2)**8 c tec7=q7(r,s,XL)-q7(r,s,-XL) uec7=z7(r,s,XL)-z7(r,s,-XL) c upot9=1.-(R1/R2)**9 c tec8=q8(r,s,XL)-q8(r,s,-XL) uec8=z8(r,s,XL)-z8(r,s,-XL) c upot10=1.-(R1/R2)**10 c tec9=q9(r,s,XL)-q9(r,s,-XL) uec9=z9(r,s,XL)-z9(r,s,-XL) c upot11=1.-(R1/R2)**11 c VECTOR POTENTIAL vector(k,i)=b00/r*R2**3*(tec1*upot3/6. + +R2**2*(-tec3*upot5/80. + +R2**2*(tec5*upot7/2688. + +R2**2*(-tec7*upot9/165888. + +R2**2*tec9*upot11/8110080.)))) c there is factor of 1/r in the functions tec* which is explicitly taken out c RADIAL MAGNETIC FIELD fldr(k,i)=b00/r*R2**3*(tec2*upot3/6. + +R2**2*(-tec4*upot5/80. + +R2**2*(tec6*upot7/2688. + +R2**2*(-tec8*upot9/165888.)))) c LONG MAGNETIC FIELD flds(k,i)=-b00*R2**3*(upot3*(uec1)/6. + +R2**2*(-upot5*(uec3)/80. + +R2**2*(upot7*(uec5)/2688. + +R2**2*(-upot9*(uec7)/165888. + +R2**2*upot11*(uec9)/8110080.)))) c if(i.eq.imiddle)then write(40,*)rad(k),vector(k,i) write(41,*)rad(k),fldr(k,i) write(42,*)rad(k),flds(k,i) endif s=s+ds end do if(k.eq.1)then do i=1,n_long_points+1 write(30,*)xlong(i),vector(k,i) write(31,*)xlong(i),fldr(k,i) write(32,*)xlong(i),flds(k,i) end do else if(k.eq.kmiddle)then do i=1,n_long_points+1 write(30,*)xlong(i),vector(k,i) write(31,*)xlong(i),fldr(k,i) write(32,*)xlong(i),flds(k,i) end do else if(k.eq.n_radial_points+1)then do i=1,n_long_points+1 write(30,*)xlong(i),vector(k,i) write(31,*)xlong(i),fldr(k,i) write(32,*)xlong(i),flds(k,i) end do endif endif ! end testing radial position r=r+dr end do write(29)n_radial_points,n_long_points write(29)vector write(28)n_radial_points,n_long_points write(28)fldr write(27)n_radial_points,n_long_points write(27)flds write(6,*)'pot_sheet.dat contains:A(i,j)', + n_radial_points+1,n_long_points+1 write(6,*)'br_sheet.dat contains:A(i,j)', + n_radial_points+1,n_long_points+1 write(6,*)'bs_sheet.dat contains:A(i,j)', + n_radial_points+1,n_long_points+1 write(6,*)'fort.30 contains: s, vector(100,s)' write(6,*)'fort.40 contains: r, vector(r,250)' write(6,*)'fort.31 contains: s, br(100,s)' write(6,*)'fort.41 contains: r, br(r,250)' write(6,*)'fort.32 contains: s, bs(100,s)' write(6,*)'fort.42 contains: r, bs(r,250)' c 50 format(1x,6(g12.5,1x)) stop end c FUNCTIONS c real*8 function a1(ss,R,S) ! first order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 c aux0=(S - ss)*Log(R + Sqrt(R**2 + (-S + ss)**2)) a1=aux0 end c real*8 function a2(ss,r,S) ! second order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0=-((ss - S)**2/ - ((r + Sqrt(r**2 + (ss - S)**2))* - Sqrt(r**2 + (ss - S)**2))) - - Log(r + Sqrt(r**2 + (ss - S)**2)) a2=aux0 end real*8 function a3(ss,r,S) ! 3rd order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 c aux0= -(((3*r**3 + 3*r**2*Sqrt(r**2 + (ss - S)**2) + - 2*r*(ss - S)**2 + - Sqrt(r**2 + (ss - S)**2)*(ss - S)**2)*(ss - S))/ - ((r + Sqrt(r**2 + (s - ss)**2))**2* - (r**2 + (ss - S)**2)**1.5)) a3=aux0 end c real*8 function a4(ss,r,S) ! 4th order implicit double precision (a-h,o-z) c eta is position of center and R1 radius common/parameters/XL,R1,R2 aux0= (-6*r**6 - 6*r**5*Sqrt(r**2 + (ss - S)**2) + - 3*r**4*(ss - S)**2 + - 6*r**3*Sqrt(r**2 + (ss - S)**2)*(ss - S)**2 + - 7*r**2*(ss - S)**4 + - 3*r*Sqrt(r**2 + (ss - S)**2)*(ss - S)**4 + (ss - S)**6 - )/((r + Sqrt(r**2 + (ss - S)**2))**3* - (r**2 + (ss - S)**2)**2.5) a4=aux0 end c real*8 function a5(ss,r,S) ! 5th order implicit double precision (a-h,o-z) c eta is position of center and R1 radius common/parameters/XL,R1,R2 aux0=((90*r**7 + 90*r**6*Sqrt(r**2 + (ss - S)**2) + - 55*r**5*(ss - S)**2 + - 10*r**4*Sqrt(r**2 + (ss - S)**2)*(ss - S)**2 - - 28*r**3*(ss - S)**4 - - 22*r**2*Sqrt(r**2 + (ss - S)**2)*(ss - S)**4 - - 8*r*(ss - S)**6 - - 2*Sqrt(r**2 + (ss - S)**2)*(ss - S)**6)*(ss - S))/ - ((r + Sqrt(r**2 + (ss - S)**2))**4*(r**2 + (ss - S)**2)**3.5) a5=aux0 end c real*8 function a6(ss,r,S) ! 6th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0=(3*(60*r**10 + 60*r**9*Sqrt(r**2 + (ss - S)**2) - _ 355*r**8*(ss - S)**2 - - 385*r**7*Sqrt(r**2 + (ss - S)**2)*(ss - S)**2 - _ 440*r**6*(ss - ss)**4 - - 240*r**5*Sqrt(r**2 + (ss - S)**2)*(ss - S)**4 - _ 28*r**4*(ss - S)**6 + - 40*r**3*Sqrt(r**2 + (ss - S)**2)*(ss - S)**6 + _ 34*r**2*(ss - S)**8 + - 10*r*Sqrt(r**2 + (ss - S)**2)*(ss - S)**8 + _ 2*(ss - S)**10))/ - ((r + Sqrt(r**2 + (ss - S)**2))**5*(r**2 + (ss - S)**2)**4.5) a6=aux0 end c real*8 function a7(s,r,ss) ! 7th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (-3*(2800*r**11 + 2800*r**10*Sqrt(r**2 + (s - ss)**2) - - 2450*r**9*(s - ss)**2 - - 3850*r**8*Sqrt(r**2 + (s - ss)**2)*(s - ss)**2 - - 7707*r**7*(s - ss)**4 - - 5432*r**6*Sqrt(r**2 + (s - ss)**2)*(s - ss)**4 - - 2556*r**5*(s - ss)**6 - - 496*r**4*Sqrt(r**2 + (s - ss)**2)*(s - ss)**6 + - 264*r**3*(s - ss)**8 + - 184*r**2*Sqrt(r**2 + (s - ss)**2)*(s - ss)**8 + - 48*r*(s - ss)**10 + - 8*Sqrt(r**2 + (s - ss)**2)*(s - ss)**10)*(s - ss))/ - ((r + Sqrt(r**2 + (s - ss)**2))**6*(r**2 + (s - ss)**2)**5.5) a7=aux0 end c real*8 function a8(s,r,ss) ! 8th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0=(-15*(1120*r**14 + 1120*r**13*Sqrt(r**2 + (s - ss)**2) - - 17220*r**12*(s - ss)**2 - - 17780*r**11*Sqrt(r**2 + (s - ss)**2)*(s - ss)**2 - - 9919*r**10*(s - ss)**4 - - 889*r**9*Sqrt(r**2 + (s - ss)**2)*(s - ss)**4 + - 21504*r**8*(s - ss)**6 + - 19656*r**7*Sqrt(r**2 + (s - ss)**2)*(s - ss)**6 + - 15000*r**6*(s - ss)**8 + - 6216*r**5*Sqrt(r**2 + (s - ss)**2)*(s - ss)**8 + - 984*r**4*(s - ss)**10 - - 336*r**3*Sqrt(r**2 + (s - ss)**2)*(s - ss)**10 - - 248*r**2*(s - ss)**12 - - 56*r*Sqrt(r**2 + (s - ss)**2)*(s - ss)**12 - - 8*(s - ss)**14))/ - ((r + Sqrt(r**2 + (s - ss)**2))**7*(r**2 + (s - ss)**2)**6.5) a8=aux0 end c real*8 function a9(s,r,ss) ! 8th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0=(15*(-(((r + Sqrt(r**2 + (s - ss)**2))* - (-34440*r**12 + (1120*r**13)/Sqrt(r**2 + (s - ss)**2) - - 35560*r**11*Sqrt(r**2 + (s - ss)**2) - - 39676*r**10*(s - ss)**2 - - (17780*r**11*(s - ss)**2)/Sqrt(r**2 + (s - ss)**2) - - 3556*r**9*Sqrt(r**2 + (s - ss)**2)*(s - ss)**2 + - 129024*r**8*(s - ss)**4 - - (889*r**9*(s - ss)**4)/Sqrt(r**2 + (s - ss)**2) + - 117936*r**7*Sqrt(r**2 + (s - ss)**2)*(s - ss)**4 + - 120000*r**6*(s - ss)**6 + - (19656*r**7*(s - ss)**6)/Sqrt(r**2 + (s - ss)**2) + - 49728*r**5*Sqrt(r**2 + (s - ss)**2)*(s - ss)**6 + - 9840*r**4*(s - ss)**8 + - (6216*r**5*(s - ss)**8)/Sqrt(r**2 + (s - ss)**2) - - 3360*r**3*Sqrt(r**2 + (s - ss)**2)*(s - ss)**8 - - 2976*r**2*(s - ss)**10 - - (336*r**3*(s - ss)**10)/Sqrt(r**2 + (s - ss)**2) - - 672*r*Sqrt(r**2 + (s - ss)**2)*(s - ss)**10 - - 112*(s - ss)**12 - - (56*r*(s - ss)**12)/Sqrt(r**2 + (s - ss)**2)))/ - (r**2 + (s - ss)**2)**6.5) + - (13*(r + Sqrt(r**2 + (s - ss)**2))* - (1120*r**14 + 1120*r**13*Sqrt(r**2 + (s - ss)**2) - - 17220*r**12*(s - ss)**2 - - 17780*r**11*Sqrt(r**2 + (s - ss)**2)*(s - ss)**2 - - 9919*r**10*(s - ss)**4 - - 889*r**9*Sqrt(r**2 + (s - ss)**2)*(s - ss)**4 + - 21504*r**8*(s - ss)**6 + - 19656*r**7*Sqrt(r**2 + (s - ss)**2)*(s - ss)**6 + - 15000*r**6*(s - ss)**8 + - 6216*r**5*Sqrt(r**2 + (s - ss)**2)*(s - ss)**8 + - 984*r**4*(s - ss)**10 - - 336*r**3*Sqrt(r**2 + (s - ss)**2)*(s - ss)**10 - - 248*r**2*(s - ss)**12 - - 56*r*Sqrt(r**2 + (s - ss)**2)*(s - ss)**12 - - 8*(s - ss)**14))/(r**2 + (s - ss)**2)**7.5 + - (7*(1120*r**14 + 1120*r**13*Sqrt(r**2 + (s - ss)**2) - - 17220*r**12*(s - ss)**2 - - 17780*r**11*Sqrt(r**2 + (s - ss)**2)*(s - ss)**2 - - 9919*r**10*(s - ss)**4 - - 889*r**9*Sqrt(r**2 + (s - ss)**2)*(s - ss)**4 + - 21504*r**8*(s - ss)**6 + - 19656*r**7*Sqrt(r**2 + (s - ss)**2)*(s - ss)**6 + - 15000*r**6*(s - ss)**8 + - 6216*r**5*Sqrt(r**2 + (s - ss)**2)*(s - ss)**8 + - 984*r**4*(s - ss)**10 - - 336*r**3*Sqrt(r**2 + (s - ss)**2)*(s - ss)**10 - - 248*r**2*(s - ss)**12 - - 56*r*Sqrt(r**2 + (s - ss)**2)*(s - ss)**12 - - 8*(s - ss)**14))/(r**2 + (s - ss)**2)**7)* - (s - ss))/(r + Sqrt(r**2 + (s - ss)**2))**8 a9=aux0 end c real*8 function q1(r,ss,S) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0=(-S + ss)/Sqrt(r**2 + (-S + ss)**2) q1=aux0 end c real*8 function q2(r,ss,S) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= -(((S - ss)*(-S + ss))/ - ((r**2 + (-S + ss)**2)**1.5)) - - 1/(Sqrt(r**2 + (-S + ss)**2)) q2=aux0 end c real*8 function q3(r,ss,S) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (3*(S - ss)*(-S + ss)**2)/ - ((r**2 + (-S + ss)**2)**2.5) - - (S - ss)/((r**2 + (-S + ss)**2)**1.5) + - (2*(-S + ss))/((r**2 + (-S + ss)**2)**1.5) q3=aux0 end c real*8 function q4(r,ss,S) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (-15*(S - ss)*(-S + ss)**3)/ - (r*(r**2 + (-S + ss)**2)**3.5) + - (9*(S - ss)*(-S + ss))/(r*(r**2 + (-S + ss)**2)**2.5) - - (9*(-S + ss)**2)/(r*(r**2 + (-S + ss)**2)**2.5) + - 3/(r*(r**2 + (-S + ss)**2)**1.5) q4=aux0 end c real*8 function q5(r,ss,S) ! 5th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (105*(S - ss)*(-S + ss)**4)/ - (r*(r**2 + (-S + ss)**2)**4.5) - - (90*(S - ss)*(-S + ss)**2)/ - (r*(r**2 + (-S + ss)**2)**3.5) + - (60*(-S + ss)**3)/(r*(r**2 + (-S + ss)**2)**3.5) + - (9*(S - ss))/(r*(r**2 + (-S + ss)**2)**2.5) - - (36*(-S + ss))/(r*(r**2 + (-S + ss)**2)**2.5) q5=aux0 end c real*8 function q6(r,ss,S) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (-945*(S - ss)*(-S + ss)**5)/ - (r*(r**2 + (-S + ss)**2)**5.5) + - (1050*(S - ss)*(-S + ss)**3)/ - (r*(r**2 + (-S + ss)**2)**4.5) - - (525*(-S + ss)**4)/(r*(r**2 + (-S + ss)**2)**4.5) - - (225*(S - ss)*(-S + ss))/ - (r*(r**2 + (-S + ss)**2)**3.5) + - (450*(-S + ss)**2)/(r*(r**2 + (-S + ss)**2)**3.5) - - 45/(r*(r**2 + (-S + ss)**2)**2.5) q6=aux0 end c real*8 function q7(r,ss,S) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (10395*(S - ss)*(-S + ss)**6)/ - (r*(r**2 + (-S + ss)**2)**6.5) - - (14175*(S - ss)*(-S + ss)**4)/ - (r*(r**2 + (-S + ss)**2)**5.5) + - (5670*(-S + ss)**5)/(r*(r**2 + (-S + ss)**2)**5.5) + - (4725*(S - ss)*(-S + ss)**2)/ - (r*(r**2 + (-S + ss)**2)**4.5) - - (6300*(-S + ss)**3)/(r*(r**2 + (-S + ss)**2)**4.5) - - (225*(S - ss))/(r*(r**2 + (-S + ss)**2)**3.5) + - (1350*(-S + ss))/(r*(r**2 + (-S + ss)**2)**3.5) q7=aux0 end c real*8 function q8(r,ss,S) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (-135135*(S - ss)*(-S + ss)**7)/ - (r*(r**2 + (-S + ss)**2)**7.5) + - (218295*(S - ss)*(-S + ss)**5)/ - (r*(r**2 + (-S + ss)**2)**6.5) - - (72765*(-S + ss)**6)/(r*(r**2 + (-S + ss)**2)**6.5) - - (99225*(S - ss)*(-S + ss)**3)/ - (r*(r**2 + (-S + ss)**2)**5.5) + - (99225*(-S + ss)**4)/(r*(r**2 + (-S + ss)**2)**5.5) + - (11025*(S - ss)*(-S + ss))/ - (r*(r**2 + (-S + ss)**2)**4.5) - - (33075*(-S + ss)**2)/(r*(r**2 + (-S + ss)**2)**4.5) + - 1575/(r*(r**2 + (-S + ss)**2)**3.5) q8=aux0 end c real*8 function q9(r,ss,s) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (2027025*(S - ss)*(-S + ss)**8)/ - (r*(r**2 + (-S + ss)**2)**8.5) - - (3783780*(S - ss)*(-S + ss)**6)/ - (r*(r**2 + (-S + ss)**2)**7.5) + - (1081080*(-S + ss)**7)/(r*(r**2 + (-S + ss)**2)**7.5) + - (2182950*(S - ss)*(-S + ss)**4)/ - (r*(r**2 + (-S + ss)**2)**6.5) - - (1746360*(-S + ss)**5)/(r*(r**2 + (-S + ss)**2)**6.5) - - (396900*(S - ss)*(-S + ss)**2)/ - (r*(r**2 + (-S + ss)**2)**5.5) + - (793800*(-S + ss)**3)/(r*(r**2 + (-S + ss)**2)**5.5) + - (11025*(S - ss))/(r*(r**2 + (-S + ss)**2)**4.5) - - (88200*(-S + ss))/(r*(r**2 + (-S + ss)**2)**4.5) q9=aux0 end c c real*8 function z1(r,ss,s) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0=(-S + ss)/(r**2 + (-S + ss)**2)**1.5 z1=aux0 end c real*8 function z2(r,ss,s) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (-3*(-S + ss)**2)/(r**2 + (-S + ss)**2)**2.5 + - (r**2 + (-S + ss)**2) z2=aux0 end c real*8 function z3(r,ss,s) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (15*(-S + ss)**3)/(r**2 + (-S + ss)**2)**3.5 - - (9*(-S + ss))/(r**2 + (-S + ss)**2)**2.5 z3=aux0 end c real*8 function z4(r,ss,s) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (-105*(-S + ss)**4)/(r**2 + (-S + ss)**2)**4.5 + - (90*(-S + ss)**2)/(r**2 + (-S + ss)**2)**3.5 - - 9/(r**2 + (-S + ss)**2)**2.5 z4=aux0 end c real*8 function z5(r,ss,s) ! 5th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (945*(-S + ss)**5)/(r**2 + (-S + ss)**2)**5.5 - - (1050*(-S + ss)**3)/(r**2 + (-S + ss)**2)**4.5 + - (225*(-S + ss))/(r**2 + (-S + ss)**2)**3.5 z5=aux0 end c real*8 function z6(r,ss,s) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (-10395*(-S + ss)**6)/(r**2 + (-S + ss)**2)**6.5 + - (14175*(-S + ss)**4)/(r**2 + (-S + ss)**2)**5.5 - - (4725*(-S + ss)**2)/(r**2 + (-S + ss)**2)**4.5 + - 225/(r**2 + (-S + ss)**2)**3.5 z6=aux0 end c real*8 function z7(r,ss,s) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (135135*(-S + ss)**7)/(r**2 + (-S + ss)**2)**7.5 - - (218295*(-S + ss)**5)/(r**2 + (-S + ss)**2)**6.5 + - (99225*(-S + ss)**3)/(r**2 + (-S + ss)**2)**5.5 - - (11025*(-S + ss))/(r**2 + (-S + ss)**2)**4.5 z7=aux0 end c real*8 function z8(r,ss,s) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (-2027025*(-S + ss)**8)/(r**2 + (-S + ss)**2)**8.5 + - (3783780*(-S + ss)**6)/(r**2 + (-S + ss)**2)**7.5 - - (2182950*(-S + ss)**4)/(r**2 + (-S + ss)**2)**6.5 + - (396900*(-S + ss)**2)/(r**2 + (-S + ss)**2)**5.5 - - 11025/(r**2 + (-S + ss)**2)**4.5 z8=aux0 end c real*8 function z9(r,ss,s) ! 1th order implicit double precision (a-h,o-z) common/parameters/XL,R1,R2 aux0= (34459425*(-S + ss)**9)/(r**2 + (-S + ss)**2)**9.5 - - (72972900*(-S + ss)**7)/(r**2 + (-S + ss)**2)**8.5 + - (51081030*(-S + ss)**5)/(r**2 + (-S + ss)**2)**7.5 - - (13097700*(-S + ss)**3)/(r**2 + (-S + ss)**2)**6.5 + - (893025*(-S + ss))/(r**2 + (-S + ss)**2)**5.5 z9=aux0 end c