求助！！！防波堤设计计算，波浪折射算角度的程序

l3319799 · 发表于 2014-4-30 16:53:02

本帖最后由 l3319799 于 2014-4-30 17:08 编辑

这是我是程序，大体思路是把海图手工打上网格做成矩阵，再经行插值各点水深，求角度
现在问题就是每次运行最后都有个-COS：或者-MOD：domain error，求解释，无错误，无警告，头要炸了

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

program main
implicit none
integer::r,i,j,m,n,k,l,h
real::T0,s,g,pi,SL,dalpha1,dalpha2,depth,L1,L0,c0,c1,p,q,beta,gamma,Lsy,theta1,theta2,theta0
!!!!!导入水深矩阵
real,dimension(0:16,0:8)::d
real,dimension(0:1000)::x,y
real,dimension(1:1000)::ZZJG
open(1,file='matrix5.dat')
do j=0,8
    do i=0,16
  read(1,*),d(i,j)
  write (*,*) d(i,j)
enddo
enddo
close(1)
!!!!
!!!起点各种参数
print *,"please enter:T0,start position m=，n=,入射角theta0="
read *,T0,m,n,s
g=9.81
pi=3.1415926
SL=200.0
L0=g*(T0**2)/(2*pi)
c0=L0/T0
x(0)=m*400.0
y(0)=n*600.0  
theta0=s
!!!从起点算至工程点
do r=0,1000
    theta0=theta0
    x(r+1)=x(r)+cosd(theta0)*SL                                                        !!!!!!!!!!!!!!!theta是角度
    y(r+1)=y(r)+sind(theta0)*SL
    p=x(r+1)/400.0
    q=y(r+1)/600.0
!!!!!判断区间
k=0
l=0
    do i=0,16
      k=k
   beta=p-k
    if (beta>m) then
       k=i+1
    else 
    endif
    enddo
    do j=0,8
      l=l
   gamma=q-l
     if (gamma>n) then
       l=j+1
     else
  endif
    enddo
!!!!!
         dalpha1=abs((d(k-1,l)-d(k-1,l-1)))*(q-(l-1))+min(d(k-1,l-1),d(k-1,l))
   dalpha2=abs((d(k,l)-d(k,l-1)))*(q-(l-1))+min(d(k,l-1),d(k,l))
   depth=abs(dalpha1-dalpha2)*(p-(k-1))+min(dalpha1,dalpha2)
   if (depth<=1.2) then 
   go to 33
      else if (depth>1.2) then
!!!!!水深插值
!!!!!
                         L1=Lsy(depth)
                   c1=L1/T0
        !角度
                theta1=theta0+((y(r+1)-y(r))*2/(c1+c0))*((c1-c0)/(x(r+1)-x(r))-(1/tan(theta0))*((c1-c0)/(y(r+1)-y(r)))) !角度
                   ZZJG(r+1)=theta1 !角度
                theta0=theta1
       c0=c1 
   ENDIF
   
!!!!!角度计算
enddo
33 do h=1,1000
print *,"各次角度",ZZJG(h)
enddo
endprogram

function Lsy(dd) result(finals)
implicit none
integer::t
real::eps,dd,temp,finals,L0,g,T0,pi
real,dimension(0:1000)::Lr
eps=0.01
Lr(0)=L0
do t=0,1000
Lr(t+1)=g*(T0**2)*tanh(2*pi*dd/L0)/(2*pi)
temp=Lr(t+1)-Lr(t)
if ((abs(temp))>=eps) then
L0=Lr(t+1)
else if (abs(temp)<eps) then
finals=Lr(t+1)
endif
end do
endfunction