[Fortran] 纯文本查看 复制代码
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