Fortran Coder

查看: 9712|回复: 1
打印 上一主题 下一主题

[数值问题] 求解运行时的问题。daltA1没算进去

[复制链接]

17

帖子

5

主题

0

精华

入门

F 币
86 元
贡献
52 点
跳转到指定楼层
楼主
发表于 2014-5-13 11:28:35 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
我的程序没有出现错误,运行的时候有两个问题:
1.不知道为啥daltA1没算进去
2.A3和daltA1循环的时候位置对调,结果就不一样
谢谢大家关注,以下是我的程序
[Fortran] 纯文本查看 复制代码
program tu5
 implicit none
 real*8 A1,A2,A4
 real*8 Z0,Zr,Zi
 real*8 A,B,C,D
 real*8 E
 real*8 F
 real*8 ref,Rf
 real*8 C1,D1,E1,F1,refs,t
 real,parameter::pi=3.14159
 real*8::A3=0.0
 integer i
 real*8::daltA1=0.0
 integer k
 character(len=80)::filename="data.txt"
 integer,parameter::fileid=10

 !!!!!!!!!!
 A2=0.25
 A4=0.425
 A1=0.625
 Z0=50
 Zr=2
 Zi=-3
 !!!!!!!!!!************************************************
 open(10,file="data.txt")
    do k=1,100
       daltA1=daltA1+0.001
       do i=1,int(0.2/0.001)
          A3=A3+0.001

    call Get_Rf(A3,daltA1,A1,A2,A4,Zi,Zr,Z0,t,E,F,E1,F1,ref,refs,Rf)
   write(*,*)t
 write(*,"('A3:'F8.5,'Rf:'F8.5,'daltA1:'F8.5)")A3,Rf,daltA1
 write(10,"('A3:'F8.5,'Rf:'F8.5,'daltA1:'F8.5)")A3,Rf,daltA1

 end do
 end do
 close(fileid)
 stop
 end
 !!!!!!!!***************************************************

 !!!!!!!!!!!!!***********************************************求ref(A1)*****************
subroutine A3_to_A(A3,A2,A4,A)
 implicit none
 real*8 A2,A3,A4,A
 real,parameter::pi=3.14159

 A=cos(2.0*pi*A2)*cos(2.0*pi*A4)-sin(2.0*pi*A2)*(-cos(2.0*pi*A4)/tan(2.0*pi*A3)+sin(2.0*pi*A4))              !!A的函数式
!write(*,*)A
 return
 end

 subroutine A3_to_B(A3,A2,A4,Z0,B)
 implicit none
 real*8 A2,A4,B,Z0,A3
 real,parameter::pi=3.14159


 B=Z0*(cos(2.0*pi*A2)*sin(2.0*pi*A4)+sin(2.0*pi*A2)*(sin(2.0*pi*A4)/tan(2.0*pi*A3)+cos(2.0*pi*A4)))            !!B的函数式
!write(*,*)B
 return
 end
 subroutine A3_to_C(A3,A1,A2,A4,Z0,C)
 implicit none
 real*8 A1,A2,A3,A4,C,Z0
 real,parameter::pi=3.14159


 C=1.0/Z0*(cos(2.0*pi*A4)*((-cos(2.0*pi*A2))/tan(2.0*pi*A1)+sin(2.0*pi*A2))+(sin(2.0*pi*A2)/tan(2.0*pi*A1)+cos(2.0*pi*A2))*((-cos(2.0*pi*A4))/tan(2.0*pi*A3)+sin(2.0*pi*A4)))
 !write(*,*)C
 return
 end

 subroutine A3_to_D(A3,A1,A2,A4,D)
 implicit none
 real*8 A1,A2,A3,A4,D
 real,parameter::pi=3.14159

 D=-sin(2.0*pi*A4)*(-cos(2.0*pi*A2)/tan(2.0*pi*A1)+sin(2.0*pi*A2))+(sin(2.0*pi*A2)/tan(2.0*pi*A1)+cos(2.0*pi*A2))*(sin(2.0*pi*A4)/tan(2.0*pi*A3)+cos(2.0*pi*A4))
 !write(*,*)D
 return
 end
 !!!!!!!!!!!***************************************E,F,ref*******
 subroutine Get_E(Z0,Zi,Zr,A3,A1,A2,A4,E)
 implicit none
 real*8 A,B,C,D,Z0,Zi,Zr,A3,E,A1,A2,A4

 call A3_to_A(A3,A2,A4,A)
 call A3_to_B(A3,A2,A4,Z0,B)
 call A3_to_C(A3,A1,A2,A4,Z0,C)
 call A3_to_D(A3,A1,A2,A4,D)


 E=(A*Zr*(D-C*Zi)+C*Zr*(B+A*Zi))/(((D-C*Zi)**2+Zr**2*C**2)*Z0)
 !write(*,*)E
 return
 end

 subroutine Get_F(A3,Z0,Zi,Zr,A1,A2,A4,F)
 implicit none
 real*8 A,B,C,D,Zi,Zr,Z0,F,A3,A1,A2,A4

 call A3_to_A(A3,A2,A4,A)
 call A3_to_B(A3,A2,A4,Z0,B)
 call A3_to_C(A3,A1,A2,A4,Z0,C)
 call A3_to_D(A3,A1,A2,A4,D)

 F=((B+A*Zi)*(D-C*Zi)-A*C*Zr**2)/(((D-C*Zi)**2+Zr**2*C**2)*Z0)
 !write(*,*)F
 return
 end

 subroutine Get_ref(A3,A1,A2,A4,Z0,Zi,Zr,E,F,ref)
 implicit none
 real*8 E,F,A3,ref,Z0,Zi,Zr,A1,A2,A4
 call Get_E(Z0,Zi,Zr,A3,A1,A2,A4,E)
 call Get_F(A3,Z0,Zi,Zr,A1,A2,A4,F)

 ref=((E**2-1+F**2)**2+4*F**2)/(((E+1)**2+F**2)**2)
 !write(*,*)ref
 return
 end
 !!!!!!!!!!!!!!!!!!**********************************************************
 !!!!!!!!!!!!!!!!!!!*******************求ref(A1+daltA1)********************
subroutine daltA1_t0_t(daltA1,A1,t)
 implicit none
 real*8 daltA1,A1,t

 t=daltA1+A1

 !write(*,*)t
 return
 end

 subroutine Get_C1(A3,daltA1,A1,A2,A4,Z0,t,C1)
 implicit none
 real*8 A1,A2,A3,A4,C1,Z0,daltA1,t
 real,parameter::pi=3.14159

 call daltA1_t0_t(daltA1,A1,t)

 C1=1.0/Z0*(cos(2.0*pi*A4)*((-cos(2.0*pi*A2))/tan(2.0*pi*t)+sin(2.0*pi*A2))+(sin(2.0*pi*A2)/tan(2.0*pi*t)+cos(2.0*pi*A2))*((-cos(2.0*pi*A4))/tan(2.0*pi*A3)+sin(2.0*pi*A4)))
 !write(*,*)C1
 return
 end

 subroutine Get_D1(A3,daltA1,A1,A2,A4,t,D1)
 implicit none
 real*8 A1,A2,A3,A4,daltA1,t,D1
 real,parameter::pi=3.14159

 call daltA1_t0_t(daltA1,A1,t)

 D1=-sin(2.0*pi*A4)*(-cos(2.0*pi*A2)/tan(2.0*pi*t)+sin(2.0*pi*A2))+(sin(2.0*pi*A2)/tan(2.0*pi*t)+cos(2.0*pi*A2))*(sin(2.0*pi*A4)/tan(2.0*pi*A3)+cos(2.0*pi*A4))
 !write(*,*)D1
 return
 end
 !!!!!!!!!!******************E1,F1,refs**********************
 subroutine Get_E1(A3,daltA1,A1,A2,A4,Z0,Zi,Zr,t,E1)
 implicit none
 real*8 A,B,C1,D1,Z0,Zi,Zr,A3,E1,daltA1,t,A1,A2,A4

 call A3_to_A(A3,A2,A4,A)
 call A3_to_B(A3,A2,A4,Z0,B)
 call Get_C1(A3,daltA1,A1,A2,A4,Z0,t,C1)
 call Get_D1(A3,daltA1,A1,A2,A4,t,D1)


 E1=(A*Zr*(D1-C1*Zi)+C1*Zr*(B+A*Zi))/(((D1-C1*Zi)**2+Zr**2*C1**2)*Z0)

 return
 end
 subroutine Get_F1(A3,daltA1,A1,A2,A4,Zi,Zr,Z0,t,F1)
 implicit none
 real*8 A,B,C1,D1,Zi,Zr,Z0,F1,A3,daltA1,t,A1,A2,A4

 call A3_to_A(A3,A2,A4,A)
 call A3_to_B(A3,A2,A4,Z0,B)
 call Get_C1(A3,daltA1,A1,A2,A4,Z0,t,C1)
 call Get_D1(A3,daltA1,A1,A2,A4,t,D1)


 F1=((B+A*Zi)*(D1-C1*Zi)-A*C1*Zr**2)/(((D1-C1*Zi)**2+Zr**2*C1**2)*Z0)

 return
 end

 subroutine Get_refs(A3,A1,A2,A4,Zi,Zr,Z0,E1,F1,t,refs)
 implicit none
 real*8 E1,F1,A3,refs,daltA1,A1,A2,A4,Z0,Zi,Zr,t
 call Get_E1(A3,daltA1,A1,A2,A4,Z0,Zi,Zr,t,E1)
 call Get_F1(A3,daltA1,A1,A2,A4,Zi,Zr,Z0,t,F1)

 refs=((E1**2-1+F1**2)**2+4*F1**2)/(((E1+1)**2+F1**2)**2)
 !write(*,*)refs
 return
 end
 !!!!!!!!!!!!!!!!!!!!!!!!*********************************************
 subroutine Get_Rf(A3,daltA1,A1,A2,A4,Zi,Zr,Z0,t,E,F,E1,F1,ref,refs,Rf)
 implicit none
 real*8 ref,refs,Rf
 real*8 A3,daltA1,E,F,E1,F1,A1,A2,A4,Z0,Zi,Zr,t
 call Get_ref(A3,A1,A2,A4,Z0,Zi,Zr,E,F,ref)
 call Get_refs(A3,A1,A2,A4,Zi,Zr,Z0,E1,F1,t,refs)

 Rf=refs/ref

 return
 end

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

736

帖子

4

主题

0

精华

大师

农村外出务工人员

F 币
700 元
贡献
359 点

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

沙发
发表于 2014-5-13 11:33:05 | 只看该作者
让别人熟悉你的代码是一件困难的事情,何况没有任何注释。

第一
do k=1,100
       daltA1=daltA1+0.001
       do i=1,int(0.2/0.001)
          A3=A3+0.001
这里的 daltA1 和 A3 位于不同的循环嵌套内,更换位置,自然会不同,这很正常。

第二.
什么叫 daltA1 没有算进去?我不知道你希望怎么算。
关于计算结果不符合预期,别人往往帮不了你,因为别人不知道你的“预期”。所以,学会 Debug 调试是关键。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

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

GMT+8, 2024-12-23 13:44

Powered by Tencent X3.4

© 2013-2024 Tencent

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