       program HCCCl
       real*8 II1,II2,II3,Navogadro,m1,m2,m3,m4,m5
       data hplanck/6.62606876d-34/
       data Navogadro/6.02214199d23/
       data pi/3.14159265358979311599/
       data m1,m2,m3,m4,m5/ 1.007825032,2.014101778,12.d0,
     X      34.96885271,36.96590260/
       data B1,B2,B3/5572.3,5084.2,5684.2/
       conv=hplanck/(8*pi*pi*1e6)*1e20*Navogadro*1000
       II1=conv/B1
       II2=conv/B2
       II3=conv/B3
       print *,"Enter delta"
       read *,delta
       print *,"Enter epsilon"
       read *,epsilon
       print *,"Enter guess values for r1,r2,r3"
       read *,r1,r2,r3

       c1=r1
       c2=r2
       c3=r3
10     continue
         a1=ff(r1,r2,r3,m1,m3,m5)
         a11=ff(r1+delta,r2,r3,m1,m3,m5)-a1
         a12=ff(r1,r2+delta,r3,m1,m3,m5)-a1
         a13=ff(r1,r2,r3+delta,m1,m3,m5)-a1
         a2=ff(r1,r2,r3,m2,m3,m5)
         a21=ff(r1+delta,r2,r3,m2,m3,m5)-a2
         a22=ff(r1,r2+delta,r3,m2,m3,m5)-a2
         a23=ff(r1,r2,r3+delta,m2,m3,m5)-a2
         a3=ff(r1,r2,r3,m1,m3,m4)
         a31=ff(r1+delta,r2,r3,m1,m3,m4)-a3
         a32=ff(r1,r2+delta,r3,m1,m3,m4)-a3
         a33=ff(r1,r2,r3+delta,m1,m3,m4)-a3
         b1=-(a1-II1)*delta
         b2=-(a2-II2)*delta
         b3=-(a3-II3)*delta
         dd=a11*a22*a33+a12*a23*a31+a13*a21*a32
     X     -a31*a22*a13-a32*a23*a11-a33*a21*a12
         d1=b1*a22*a33+a12*a23*b3+a13*b2*a32
     X     -b3*a22*a13-a32*a23*b1-a33*b2*a12
         d2=a11*b2*a33+b1*a23*a31+a13*a21*b3
     X     -a31*b2*a13-b3*a23*a11-a33*a21*b1
         d3=a11*a22*b3+a12*b2*a31+b1*a21*a32
     X     -a31*a22*b1-a32*b2*a11-b3*a21*a12
         c1=d1/dd
         c2=d2/dd
         c3=d3/dd
         r1=r1+c1
         r2=r2+c2
         r3=r3+c3
         print *,r1,r2,r3
       if((abs(c1/r1).gt.epsilon).or.(abs(c2/r2).gt.epsilon).or.
     X       (abs(c3/r3).gt.epsilon)) goto 10
       print *,"Final results"
       print *,r1,r2,r3
       end
       
       function ff(r1,r2,r3,m1,m2,m3)
       real*8 m1,m2,m3
       xx=(m3*r3-m1*(r1+r2)-m2*r2)/(m1+2*m2+m3)
       ff=m1*(r1+r2+xx)**2+m2*((r2+xx)**2+xx**2)+m3*(r3-xx)**2
       return
       end
