c
      program circular2
      real mu
c
 10   write(*,*)'Enter Zs'
      read(*,*)Zs
c
      mu=0.3
      slip =1.0
      R1=Zs*30.0
      R2=1.0
      m1=1000
      m2=2000      
      if(R1.ge.R2)then
        R1=R2
        m2=0
      end if
      write(*,*)
      write(*,*)'mu,slip,R1,R2 =',mu,slip,R1,R2
      write(*,*)
c
      pi=3.141592654
      c3=-1.0/(8.0*pi*(1.0-mu))
      c4=1.0-2.0*mu
c
      dR=R1/m1
      write(*,*)'m1,m2,dR1',m1,m2,dR
c
      open(9,file='u1_d',status='unknown')
      write(9,2000)m1,R1
      write(9,2000)m2,R2
 2000 format(i12,e12.4)
c
      rd=0.0
c      rr=sqrt(rd*rd+Zs*Zs)
c      rr2=rd*rd/(rr*rr)
c      u1=(2.0*C4+3.0*rr2)*rd/(rr*rr*rr)
      u1=0.0
      write(*,*)'rd,u1',rd,-pi*c3*Zs*u1
      write(9,2010)rd,-pi*c3*Zs*u1
 2010 format(2e12.4)
      sum1=0.0
      do j=1,m1
        rd=dR*j
        rr=sqrt(rd*rd+Zs*Zs)
        rr2=rd*rd/(rr*rr)
        u2=(2.0*C4+3.0*rr2)*rd/(rr*rr*rr)
        sum1=sum1+0.5*(u1+u2)
        u1=u2
c      write(*,*)'rd,u1,sum',rd,-pi*c3*Zs*u1,-pi*c3*Zs*sum1*dR
      write(9,2010)rd,-pi*c3*Zs*u1
      end do
      sum1=-pi*c3*Zs*sum1*dR
c
      if(R1.ge.R2)goto 20
      dR=(R2-R1)/m2
      write(*,*)'m2,dR',m2,dR
c
      rd=R1
      rr=sqrt(rd*rd+Zs*Zs)
      rr2=rd*rd/(rr*rr)
      u1=(2.0*C4+3.0*rr2)*rd/(rr*rr*rr)
      write(*,*)'rd,u1',rd,-pi*c3*Zs*u1
      write(9,2010)rd,-pi*c3*Zs*u1
      sum2=0.0
      do j=1,m2
        rd=R1+dR*j
        rr=sqrt(rd*rd+Zs*Zs)
        rr2=rd*rd/(rr*rr)
        u2=(2.0*C4+3.0*rr2)*rd/(rr*rr*rr)
        sum2=sum2+0.5*(u1+u2)
        u1=u2
c      write(*,*)'rd,u1,sum',rd,-pi*c3*Zs*u1,-pi*c3*Zs*sum2*dR
      write(9,2010)rd,-pi*c3*Zs*u1
      end do
      sum2=-pi*c3*Zs*sum2*dR
c
 20   write(9,2010)sum1,sum2
      write(9,2010)Zs,sum1+sum2
      close(9)
c
      write(*,*)
      write(*,*)'sum =',sum1,sum2,sum1+sum2
c
c      goto 10     
c
      stop
      end
