Fortran Coder

查看: 6849|回复: 3
打印 上一主题 下一主题

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

[复制链接]

12

帖子

3

主题

0

精华

入门

F 币
54 元
贡献
30 点
跳转到指定楼层
楼主
发表于 2014-5-13 14:29:24 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
[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

分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

712

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
607 元
贡献
311 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

沙发
发表于 2014-5-13 14:35:11 | 只看该作者
请详细描述你的问题

12

帖子

3

主题

0

精华

入门

F 币
54 元
贡献
30 点
板凳
 楼主| 发表于 2014-5-13 14:54:10 | 只看该作者
这里的x,y是广义坐标,该如何定义它们?
编译警告提示它们没有定义,但x,y是没有初值的

712

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
607 元
贡献
311 点

新人勋章爱心勋章水王勋章元老勋章热心勋章

地板
发表于 2014-5-13 15:00:10 | 只看该作者
你要学会用程序语言来描述你的问题。Fortran语法里没有广义坐标这个东西,只有real,integer这些数据类型。

你希望 X Y 是实数,那么就用 real 来定义它们。

初值的问题,也取决与你的想法。应该是什么初值,就怎么赋值。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-4-26 22:08

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表