subroutine ellipse ********************************************************************* * * * This subroutine is to solve the ellipse * * equation of the potential phi(i,j). * * * * Version 0.5 - Last Modified 5/18/04 * * * ********************************************************************* INCLUDE "fpvar.inc" dimension phi2(0:nx) ********************************************************************* * * * Give basic parameters. * * * ********************************************************************* re=2.0e-8 nit=0 dx1=1.0/(dx**2) dz1=1.0/(dz**2) c b1=-0.5/(dx1+dz1) it=1 a2=100.0 r=1.0 eps0=1.0 epc=0.0 epmax=10000.0 ********************************************************************* * * * Limit the amplitude of the potential phi(i,j). * * * ********************************************************************* phimax=15.0 phimin=-15.0 do 200, i=1,moz do 200, j=-5,nx if (phi(i,j).lt.phimin) phi(i,j)=phimin 200 if (phi(i,j).gt.phimax) phi(i,j)=phimax ********************************************************************* * * * Calculate the values of relevant quantities. * * * ********************************************************************* do 220, i=ST,nz cnx(i)=-2.0*(cosi**2)*ti(i)/(1.0+rvi(i)**2) cnz(i)=-2.0*(sini**2)*ti(i)/(1.0+rvi(i)**2) 220 dn(i)=2.0*ri2(i)*ti(i) do 222, i=ST,nz do 222, j=ST,nx fx(i,j)=rix1(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) $ -ri2(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b)) $ -rex1(i)*(uox+u1x(i,j)-eox/(rve(i)*b)) $ +re2(i)*(uoz+u1z(i,j)-eoz/(rve(i)*b)) $ -eoy*sini/((1.0+rvi(i)**2)*b) $ +eoy*sini/((1.0+rve(i)**2)*b) $ +sini*cosi*gg/((1.0+rvi(i)**2)*vin(i)) fz(i,j)=-ri2(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) $ +riz1(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b)) $ +re2(i)*(uox+u1x(i,j)-eox/(rve(i)*b)) $ -rez1(i)*(uoz+u1z(i,j)-eoz/(rve(i)*b)) $ -eoy*cosi/((1.0+rvi(i)**2)*b) $ +eoy*cosi/((1.0+rve(i)**2)*b) $ -(rvi(i)**2+sini**2)*gg/((1.0+rvi(i)**2)*vin(i)) 222 continue ********************************************************************* * * * Determine the values of phi(i,j) at the lower boundary. * * * ********************************************************************* 230 do 232, j=1,mox 232 phi2(j)=phi(2,j) do 234, j=1,mox 234 phi(1,j)=phi2(j)*r+(1.0-r)*phi(1,j) phi(1,mox+1)=phi(1,2) phi(1,0)=phi(1,mox-1) ********************************************************************* * * * Calculate phi(i,j) in the inner region (i=3:moz-2, j=1:mox) * * * ********************************************************************* eps1=0.0 do 265, i=3,moz-2 aven=0.0 do 240, j=1,mox-1 240 aven=aven+n(i,j) aven=aven/real(mox-1) do 245, j=1,mox aphix=n(i,j)*(rix1(i)/(rvi(i)*b)+rex1(i)/(rve(i)*b)) aphiz=n(i,j)*(riz1(i)/(rvi(i)*b)+rez1(i)/(rve(i)*b)) bphi=-n(i,j)*(ri2(i)/(rvi(i)*b)+re2(i)/(rve(i)*b)) ax=(rix1(i)/(rvi(i)*b)+rex1(i)/(rve(i)*b)) $ *(n(i,j+1)-n(i,j-1))/(2.0*dx*aphix) $ -(n(i+1,j)*ri2(i+1)/(rvi(i+1)*b) $ +n(i+1,j)*re2(i+1)/(rve(i+1)*b) $ -n(i-1,j)*ri2(i-1)/(rvi(i-1)*b) $ -n(i-1,j)*re2(i-1)/(rve(i-1)*b))/(2.0*dz*aphix) az=(n(i+1,j)*riz1(i+1)/(rvi(i+1)*b) $ +n(i+1,j)*rez1(i+1)/(rve(i+1)*b) $ -n(i-1,j)*riz1(i-1)/(rvi(i-1)*b) $ -n(i-1,j)*rez1(i-1)/(rve(i-1)*b))/(2.0*dz*aphix) $ -(ri2(i)/(rvi(i)*b)+re2(i)/(rve(i)*b)) $ *(n(i,j+1)-n(i,j-1))/(2.0*dx*aphix) s14=2.0*ti(i)*(cosi*xk-sini*zkr(i))*(cosi*xk-sini*zkr(i)) s14=s14/((1.0+rvi(i)*rvi(i))*aphix) s14=s14*(n(i,j)-aven)/aven s5=(n(i,j+1)*fx(i,j+1)-n(i,j-1)*fx(i,j-1))/(2.0*dx*aphix) s6=(n(i+1,j)*fz(i+1,j)-n(i-1,j)*fz(i-1,j))/(2.0*dz*aphix) sij=s14+s5+s6 dzz1=dz1*aphiz/aphix b1=-0.5/(dx1+dzz1) abij=bphi/(2.0*dx*dz*aphix) a1ij=dx1+ax/(2.0*dx) c1ij=dx1-ax/(2.0*dx) d1ij=dzz1+az/(2.0*dz) e1ij=dzz1-az/(2.0*dz) phi2(j)=b1*(sij-abij*phi(i+1,j+1)+abij*phi(i+1,j-1) $ +abij*phi(i-1,j+1)-abij*phi(i-1,j-1) $ -a1ij*phi(i,j+1)-c1ij*phi(i,j-1) $ -d1ij*phi(i+1,j)-e1ij*phi(i-1,j)) ********************************************************************* * * * Restrict the sudden change of phi2(i,j). * * * ********************************************************************* phiap=5.0*(phi(i,j+1)+phi(i,j-1)+phi(i+1,j)+phi(i-1,j)) phian=-5.0*(phi(i,j+1)+phi(i,j-1)+phi(i+1,j)+phi(i-1,j)) if (t.gt.20.0) then if (phi2(j).gt.phiap) then phi2(j)=0.25*(phi(i,j+1)+phi(i,j-1)+phi(i+1,j)+phi(i-1,j)) end if if (phi2(j).lt.phian) then phi2(j)=0.25*(phi(i,j+1)+phi(i,j-1)+phi(i+1,j)+phi(i-1,j)) end if end if 245 continue ********************************************************************* * * * Limit the amplitude of the potential phi2(j). * * * ********************************************************************* do 250, j=1,mox if (phi2(j).lt.phimin) then phi2(j)=phimin end if if (phi2(j).gt.phimax) then phi2(j)=phimax end if 250 continue ********************************************************************* * * * Put the values of phi2(j) to phi(i,j). * * * ********************************************************************* do 260, j=1,mox eps1=max(eps1,abs((phi2(j)-phi(i,j))*r)) 260 phi(i,j)=phi2(j)*r+(1.0-r)*phi(i,j) phi(i,mox+1)=phi(i,2) phi(i,0)=phi(i,mox-1) if (eps1.gt.epc) then epc=eps1 end if 265 continue ********************************************************************* * * * If 'epc' is bigger than 100, stop the calculation. * * If the calculation becomes divergent, that is, 'eps1' * * begins to increase, stop the calculation. * * * ********************************************************************* nit=nit+1 if (epc.gt.100.0) go to 280 if (eps1.lt.epmax) then epmax=eps1 else go to 280 end if ********************************************************************* * * * Determine the values of phi(2,j), phi(1,j), phi(moz-1,j), * * and phi(moz,j) for j=1,mox under suitable boundary * * conditions. These values are not calculated in the * * calculations for the inner region. * * * ********************************************************************* do 270, j=1,mox phi(2,j)=phi(4,j)-(cosi/sini)*(phi(3,j+1)-phi(3,j-1))*dz/dx phi(1,j)=phi(3,j)-(cosi/sini)*(phi(2,j+1)-phi(2,j-1))*dz/dx phi(moz-1,j)=phi(moz-3,j) $ +(cosi/sini)*(phi(moz-2,j+1)-phi(moz-2,j-1))*dz/dx phi(moz,j)=phi(moz-2,j) $ +(cosi/sini)*(phi(moz-1,j+1)-phi(moz-1,j-1))*dz/dx 270 continue phi(moz,mox+1)=phi(moz,2) phi(moz,0)=phi(moz,mox-1) ********************************************************************* * * * Calculate the relaxation factor 'r'. * * In fact, 'r' is taken to be one in the present calculations * * * ********************************************************************* c if (it.eq.1) then c a1=a2 c a2=alog10(eps1/eps0) c if (abs(a2-a1).lt.0.006) then c it=2 c if (eps1/eps0.ge.1.0) then c r=r-0.1 c else c r=2.0/(1.0+sqrt(1.0-eps1/eps0)) c end if c end if c end if c if (r.gt.1.5) r=1.5 c if (r.lt.1.0) r=1.0 r=1.0 ********************************************************************* * * * Determine the accuracy of the calculated data. * * If the absolute error 'abs(esp1-eps0)' or the relative * * error 'abs(1.0-eps0/eps1)' is smaller than the value * * required, stop calculation. * * If iteration number is more than 50, stop calculation. * * Otherwise, return to the beginning and calculate again. * * * ********************************************************************* c if (abs(eps1-eps0).lt.re) then if (abs(1.0-eps0/eps1).lt.re) go to 280 if (nit.gt.50) go to 280 c write(8,*)'eps1= ',eps1 c write(8,*)'r= ',r if (eps1.gt.re) then eps0=eps1 go to 230 end if ********************************************************************* * * * Periodic boundary condition in the horizontal direction. * * * ********************************************************************* 280 do 282, j=ST,0 do 282, i=1,moz 282 phi(i,j)=phi(i,mox+j-1) do 284, j=mox+1,nx do 284, i=1,moz 284 phi(i,j)=phi(i,j-mox+1) open (status='old',access='append',unit=8,file='trouble.dat') write(8,*)'epc= ',epc write(8,*)'nit= ',nit write(8,*)'end of ellipse' close(8) return end