[数值问题] 广义坐标的初值

program main  
implicit real(8) (a-h,o-z)

real(8) x,y,xv,yv,pi,sita,ss,pe,xx1,t1,t2,t3,a1,a2,a3
real(8) ff1,f1,ff2,f2,ff3,f3,k,xx,step,pev,faiv
real(8) ffr,fft,ffx,ffy,fai

    pi=3.14
        pe=sqrt(x**2+y**2)
    pev=(y*yv+x*xv)/pe
        faiv=(x*yv-y*xv)/pe

   ss=-(-pev/dabs(-pev))*dacos((-faiv)/(dsqrt(pev**2+faiv**2)))
    ff1=0.0
        ff2=0.0
        ff3=0.0
     t1=0.7745967
        t2=-0.7745967
        t3=0
        step=100
      dpi=(pi+ss)/step
        xx=0
        xx1=xx+dpi
                
    a1=0.5555556
        a2=0.5555556
        a3=0.8888889
    ff1=a1*cos((dpi/2)*t1+(xx1+xx)/2)**2/(1+pe*cos((dpi/2)*t1+(xx1+xx)/2))**3+a2*cos((dpi/2)*t2+(xx1+xx)/2)**2/(1+pe*cos((dpi/2)*t2+(xx1+xx)/2))**3+a3*cos((dpi/2)*t3+(xx1+xx)/2)**2/(1+pe*cos((dpi/2)*t3+(xx1+xx)/2))**3    
    f1=ff1
        do k=1,step-1
        xx=xx1
        xx1=xx1+dpi                                
  ff1=a1*cos((dpi/2)*t1+(xx1+xx)/2)**2/(1+pe*cos((dpi/2)*t1+(xx1+xx)/2))**3+a2*cos((dpi/2)*t2+(xx1+xx)/2)**2/(1+pe*cos((dpi/2)*t2+(xx1+xx)/2))**3+a3*cos((dpi/2)*t3+(xx1+xx)/2)**2/(1+pe*cos((dpi/2)*t3+(xx1+xx)/2))**3    
   f1=f1+ff1
   end do
   f1=(f1*dpi)/2
write(*,*) "f1=",f1
        
ff2=a1*dcos((dpi/2)*t1+(xx1+xx)/2)*dsin((dpi/2)*t1+(xx1+xx)/2)/(1+pe*dcos((dpi/2)*t1+(xx1+xx)/2))**3+a2*dcos((dpi/2)*t2+(xx1+xx)/2)*dsin((dpi/2)*t2+(xx1+xx)/2)/(1+pe*dcos((dpi/2)*t2+(xx1+xx)/2))**3+a3*dcos((dpi/2)*t3+(xx1+xx)/2)*dsin((dpi/2)*t3+(xx1+xx)/2)/(1+pe*dcos((dpi/2)*t3+(xx1+xx)/2))**3                                  
f2=ff2
        do k=1,step-1
        xx=xx1
        xx1=xx1+dpi
ff2=a1*dcos((dpi/2)*t1+(xx1+xx)/2)*dsin((dpi/2)*t1+(xx1+xx)/2)/(1+pe*dcos((dpi/2)*t1+(xx1+xx)/2))**3+a2*dcos((dpi/2)*t2+(xx1+xx)/2)*dsin((dpi/2)*t2+(xx1+xx)/2)/(1+pe*dcos((dpi/2)*t2+(xx1+xx)/2))**3+a3*dcos((dpi/2)*t3+(xx1+xx)/2)*dsin((dpi/2)*t3+(xx1+xx)/2)/(1+pe*dcos((dpi/2)*t3+(xx1+xx)/2))**3           
f2=f2+ff2
end do
f2=(f2*dpi)/2
!write(*,*) "f2=",f2


ff3=a1*sin((dpi/2)*t1+(xx1+xx)/2)**2/(1+pe*cos((dpi/2)*t1+(xx1+xx)/2))**3+a2*sin((dpi/2)*t2+(xx1+xx)/2)**2/(1+pe*cos((dpi/2)*t2+(xx1+xx)/2))**3+a3*sin((dpi/2)*t3+(xx1+xx)/2)**2/(1+pe*cos((dpi/2)*t3+(xx1+xx)/2))**3    

f3=ff3
        do k=1,step-1
        xx=xx1
        xx1=xx1+dpi
ff3=a1*sin((dpi/2)*t1+(xx1+xx)/2)**2/(1+pe*cos((dpi/2)*t1+(xx1+xx)/2))**3+a2*sin((dpi/2)*t2+(xx1+xx)/2)**2/(1+pe*cos((dpi/2)*t2+(xx1+xx)/2))**3+a3*sin((dpi/2)*t3+(xx1+xx)/2)**2/(1+pe*cos((dpi/2)*t3+(xx1+xx)/2))**3    
f3=f3+ff3
end do
f3=(f3*dpi)/2
!write(*,*) "f3=",f3
fai=datan(y/x)

    if (x<0) then
            fai=pi+fai
        endif

    ffr=-(f1*pev+f2*faiv)
        fft=-(f2*pev+f3*faiv)

    ffx=ffr*cos(fai)-fft*sin(fai)
        ffy=ffr*sin(fai)+fft*cos(fai)
!write(*,*) ss,y
stop
end