广义坐标的初值

liruipei · 发表于 2014-5-13 14:29:24

[Fortran] 纯文本查看 复制代码

01program main  
02implicit real(8) (a-h,o-z)
03 
04real(8) x,y,xv,yv,pi,sita,ss,pe,xx1,t1,t2,t3,a1,a2,a3
05real(8) ff1,f1,ff2,f2,ff3,f3,k,xx,step,pev,faiv
06real(8) ffr,fft,ffx,ffy,fai
07 
08    pi=3.14
09        pe=sqrt(x**2+y**2)
10    pev=(y*yv+x*xv)/pe
11        faiv=(x*yv-y*xv)/pe
12 
13   ss=-(-pev/dabs(-pev))*dacos((-faiv)/(dsqrt(pev**2+faiv**2)))
14    ff1=0.0
15        ff2=0.0
16        ff3=0.0
17     t1=0.7745967
18        t2=-0.7745967
19        t3=0
20        step=100
21      dpi=(pi+ss)/step
22        xx=0
23        xx1=xx+dpi
24                 
25    a1=0.5555556
26        a2=0.5555556
27        a3=0.8888889
28    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   
29    f1=ff1
30        do k=1,step-1
31        xx=xx1
32        xx1=xx1+dpi                                
33  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   
34   f1=f1+ff1
35   end do
36   f1=(f1*dpi)/2
37write(*,*) "f1=",f1
38         
39ff2=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                                 
40f2=ff2
41        do k=1,step-1
42        xx=xx1
43        xx1=xx1+dpi
44ff2=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          
45f2=f2+ff2
46end do
47f2=(f2*dpi)/2
48!write(*,*) "f2=",f2
49 
50 
51ff3=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   
52 
53f3=ff3
54        do k=1,step-1
55        xx=xx1
56        xx1=xx1+dpi
57ff3=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   
58f3=f3+ff3
59end do
60f3=(f3*dpi)/2
61!write(*,*) "f3=",f3
62fai=datan(y/x)
63 
64    if (x<0) then
65            fai=pi+fai
66        endif
67 
68    ffr=-(f1*pev+f2*faiv)
69        fft=-(f2*pev+f3*faiv)
70 
71    ffx=ffr*cos(fai)-fft*sin(fai)
72        ffy=ffr*sin(fai)+fft*cos(fai)
73!write(*,*) ss,y
74stop
75end