PROGRAM try
      double precision k,ka,kb,k2,k4
      logical on
      
      on = .true.
      print *,on
      stop
      
      i = 0
 50   do 100
      i = i+1
 100  while (i.le.10) 

      ONE = 1.d0
      ZERO = 0.d0
      TWO_THIRDS = 2.d0/3.d0
            
      k = 0.3
      k2 = k*k
      k4 = k2*k2
      f = rd(ZERO,1-k2,ONE)*TWO_THIRDS - rf(ZERO,1-k2,ONE)
      
      h = 0.00001
      ka = k+h
      kb = k-h
      fa = rd(ZERO,1-ka*ka,ONE)*TWO_THIRDS - rf(ZERO,1-ka*ka,ONE)

      fb = rd(ZERO,1-kb*kb,ONE)*TWO_THIRDS - rf(ZERO,1-kb*kb,ONE)
      fd = (fa-fb)/h/2.
      fd2 = (fa+fb-f-f)/h/h
      rd1 = rd(ZERO,1-k2,ONE)
      rd2 = rd(ZERO,ONE,1-k2)
      


      fdp = (-rd1*TWO_THIRDS+(2-k2)*rd2/3.)/k
      fdp2 = (rd1*(6-5*k2)-(2*k4-10*k2+6)*rd2)/3/k2/(1-k2)

      print *,'first derivative = ',fd,fdp
      print *,'second derivative = ',fd2,fdp2

      w(1) = 1.
      call xx
      
      STOP
      END

      function sum(xx)
      dimension xx(10)
      sum = 0.
      do i = 1,10
      sum = cc+xx(i)
      enddo
      return
      end
      
      subroutine xx

      print *, w(1)
      
      end
      
      FUNCTION rd(x,y,z)
      DOUBLE PRECISION rd,x,y,z,ERRTOL,TINY,BIG,C1,C2,C3,C4,C5,C6
      PARAMETER (ERRTOL=.0015d0,TINY=3D-39,BIG=1.7D38,C1=3.d0/14.d0,
     *C2=1.d0/6.d0,C3=9.d0/22.d0,C4=3.d0/26.d0,C5=.25d0*C3,C6=1.5d0*C4)
      DOUBLE PRECISION alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,
     *sqrtx,sqrty,sqrtz,sum,xt,yt,zt
      if(min(x,y).lt.0..or.min(x+y,z).lt.TINY.or.max(x,y,
     *z).gt.BIG)pause 'invalid arguments in rd'
      xt=x
      yt=y
      zt=z
      sum=0.d0
      fac=1.d0
1     continue
        sqrtx=dsqrt(xt)
        sqrty=dsqrt(yt)
        sqrtz=dsqrt(zt)
        alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz
        sum=sum+fac/(sqrtz*(zt+alamb))
        fac=.25*fac
        xt=.25*(xt+alamb)
        yt=.25*(yt+alamb)
        zt=.25*(zt+alamb)
        ave=.2*(xt+yt+3.*zt)
        delx=(ave-xt)/ave
        dely=(ave-yt)/ave
        delz=(ave-zt)/ave
      if(max(dabs(delx),dabs(dely),dabs(delz)).gt.ERRTOL)goto 1
      ea=delx*dely
      eb=delz*delz
      ec=ea-eb
      ed=ea-6.*eb
      ee=ed+ec+ec
      rd=3.*sum+fac*(1.+ed*(-C1+C5*ed-C6*delz*ee)+delz*(C2*ee+delz*(-C3*
     *ec+delz*C4*ea)))/(ave*dsqrt(ave))
      return
      END


      FUNCTION rf(x,y,z)
      DOUBLE PRECISION rf,x,y,z,ERRTOL,TINY,BIG,THIRD,C1,C2,C3,C4
      PARAMETER (ERRTOL=.0025d0,TINY=3d-39,BIG=1.7d38,THIRD=1.d0/3.d0,
     *C1=1.d0/24.d0,C2=.1d0,C3=3.d0/44.d0,C4=1.d0/14.d0)
      DOUBLE PRECISION alamb,ave,delx,dely,delz,e2,e3,sqrtx,sqrty,
     *sqrtz,xt,yt,zt
      if(min(x,y,z).lt.0..or.min(x+y,x+z,y+z).lt.TINY.or.max(x,y,
     *z).gt.BIG)pause 'invalid arguments in rf'
      xt=x
      yt=y
      zt=z
1     continue
        sqrtx=dsqrt(xt)
        sqrty=dsqrt(yt)
        sqrtz=dsqrt(zt)
        alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz
        xt=.25*(xt+alamb)
        yt=.25*(yt+alamb)
        zt=.25*(zt+alamb)
        ave=THIRD*(xt+yt+zt)
        delx=(ave-xt)/ave
        dely=(ave-yt)/ave
        delz=(ave-zt)/ave
      if(max(dabs(delx),dabs(dely),dabs(delz)).gt.ERRTOL)goto 1
      e2=delx*dely-delz**2
      e3=delx*dely*delz
      rf=(1.+(C1*e2-C2-C3*e3)*e2+C4*e3)/dsqrt(ave)
      return
      END