[updated 28.Sep.2000]
Librairie-d'application appl_rw3d > Fichier gwdyi.f |
*
subroutine gwdyi(coo,wg)
*
Calcul du champ lointain de la fonction de green
de resistance de vagues 3D
D' apres Baar & Price, 87
*
developement en serie de bessel y-i
x>0 y>0 z<0
x/d grand
*
implicit double precision (a-h)
implicit double precision (o-z)
real coo(3)
real wg(6)
dimension wgd(6)
parameter (nx =100)
parameter (nmx=2*nx+3)
parameter (qpi=12.56637061435917295)
parameter (pi=3.1415926535897)
parameter (eps=1.e-6 )
dimension byn(0:nmx),bin (0:nmx)
dimension a (0:nx) ,bn (0:nx) ,pyn (0:nx),amlt (0:nx)
dimension ax (0:nx) ,bxn (0:nx) ,pyxn(0:nx)
dimension ad (0:nx) ,bdn (0:nx)
dimension bbn (0:nx)
real bey0,bey1,bei0
na=40
naa=na
nax=na
nad=na
x=dble(coo(1))
y=dble(coo(2))
z=dble(coo(3))
d = sqrt(y*y+z*z)
ds2 = d*0.5
usd = 1/d
cb = z*usd
sb = y*usd
dcb = 2*cb
dsb = 2*sb
qpezs2= qpi*exp(z*0.5)
x2sd = x*x*usd
evaluation des ii+1/ii(x) [0:nm]
----------------------------------
if(ds2.lt.5) then
nm=int(4.*ds2+21.)
else IF(ds2.lt.10) then
nm=int(2.4*ds2+28.)
else IF(ds2.lt.20.) then
nm=int(2.*ds2+32.)
else
nm=int(1.5*ds2+42.)
endif
nm=max(nm,na)
bin(nm) =0
do 50 i=nm,1,-1
50 bin(i-1)=1./( (i+i)/ds2 + bin(i) )
by0 =dble(bey0(real(x)))
by1 =dble(bey1(real(x)))
bi0 =dble(bei0(real(ds2)))*exp(ds2)
bi1 =bin(0)*bi0
bi2 =bin(1)*bi1
pyn (0) =-1
recurence montante : an
-----------------------
a(0) = ((qpezs2+qpezs2)* by1 *bi0)
ax(0) = qpezs2 *(by2-by0)*bi0
ad(0) = qpezs2 * by1 *bi1
byn(0)=by1/by0
usbyn0=by0/by1
xsqd =x/sqrt(d)
if (xsqd.lt.2.3) then
na=int(10+6*xsqd)
else
na=int(10 *xsqd)
endif
na=25
byn(1)=(((2)/x)-1/(byn(0 )))
byn(2)=(((4)/x)-1/(byn(1)))
pyn (1)= byn(2)*byn(1)
pyxn (1)= byn(1)*byn(0 )
by2=pyxn (1)*by0
ax(0) = qpezs2*(by0-by2)*bi0
pyxn(0) =1/(byn(1)*byn(0 ))
do 200 i=1,na
ip1=i+1
i2=i+i
i4=i2+i2
byn(i2 +1)=(((i4+2)/x)-1/(byn(i2 )))
byn(i2 +2)=(((i4+4)/x)-1/(byn(i2 +1)))
pyn (ip1)=byn(i2+2)*byn(i2+1)
pyxn(ip1)=byn(i2+1)*byn(i2 )
amlt(i) =((pyn(i))*(1-(pyn(ip1)))/(1-(pyn(i))))
200 continue
amlt(0) =((pyn(0))*(1-(pyn(1)))/(1-(pyn(0))))
a (1) =(a (0)*(bin(0))*amlt(0))
ax(1) =ax(0)*bin(0)*
* (1+pyxn(1)*(pyxn(2)-2))/(1/pyxn(0)-2+pyxn(1))
ad(1) =ad(0)*amlt(0)*0.5*(bin(1)+1/bin(0))
amina=1.e6
do 210 i=1,na
ip1=i+1
a (ip1) =((a (i))*(bin(i))*amlt(i))
if ((abs(a(ip1)).gt.100*amina)
* .or.(abs(a(ip1)).lt.eps))then
naa=i
goto 220
else
amina=min(amina,abs((a(ip1))))
endif
210 continue
220 naa=min(naa,na)
aminx=1.e6
do 230 i=1,na
ip1=i+1
ax(ip1) =ax(i)*bin(i)*
* (1+pyxn(ip1)*(pyxn(i+2)-2))/(1/pyxn(i)-2+pyxn(ip1))
if ((abs(ax(ip1)).gt.100*aminx)
* .or.(abs(ax(ip1)).lt.eps))then
nax=i
goto 240
else
aminx=min(aminx,abs((ax(ip1))))
endif
230 continue
240 nax=min(nax,na)
amind=1.e6
do 250 i=1,na
ip1=i+1
ad(ip1) =ad(i)*(1+bin(ip1)*bin(i))/
* (1/bin(i-1)+bin(i))*amlt(i)
if ((abs(ad(ip1)).gt.100*amind)
* .or.(abs(ad(ip1)).lt.eps))then
nad=i
goto 260
else
amind=min(amind,abs((ad(ip1))))
endif
250 continue
260 nad=min(nad,na)
recurrence descendante : bn
---------------------------
bn (naa+1)=0
bn (naa+2)=0
bxn(nax+1)=0
bxn(nax+2)=0
bdn(nad+1)=0
bdn(nad+2)=0
bbn(naa+1)=0
bbn(naa+2)=0
nmax=max(naa,nad,nax)
do 300 i=naa,0,-1
ip1=i+1
ip2=i+2
bn (i) =(a(i)-(dcb)*(bn(ip1))-(bn(ip2)))
bbn(i) =dsb*bn(ip1)-dcb*bbn(ip1)-bbn(ip2)
300 continue
do 310 i=nax,0,-1
ip1=i+1
ip2=i+2
bxn(i) =ax(i)-dcb*bxn(ip1)-bxn(ip2)
310 continue
do 320 i=nad,0,-1
ip1=i+1
ip2=i+2
bdn(i) =ad(i)-dcb*bdn(ip1)-bdn(ip2)
320 continue
wgd(1)=((bn (0)-bn (2))*0.5)
wgd(2)=(bxn(0)-bxn(2))*0.5
wd =(bdn(0)-bdn(2))*0.5
wb =(bbn(0)-bbn(2))*0.5
wgd(5)= cb*0.5*wgd(1)+wd
wgd(6)=-d*sb*0.5*wgd(1)+wb
wgd(3)=wgd(5)*sb+cb*wgd(6)*usd
wgd(4)=wgd(5)*cb-sb*wgd(6)*usd
wg(1)=real(wgd(1))
wg(2)=real(wgd(2))
wg(3)=real(wgd(3))
wg(4)=real(wgd(4))
wg(5)=real(wgd(5))
wg(6)=real(wgd(6))
return
end
gwdyi est appelé dans