[updated 28.Sep.2000]
Librairie-d'application appl_rw3d > Fichier gwdjk.f |
*
subroutine gwdjk(coo,wg)
*
Calcul du champ lointain de la fonction de green
de resistance de vagues 3D
D' apres Baar & Price, 87
Derniere Modif. : Correction des boucles 50 et 60 (Royal Christophe).
Modif. de la declaration de bjn. (Le 06/03/1992).
Ceci corrige une division par 0 due a une erreur de
programmation. Je l'ai signalee a Eric Masson (ENSTA)
*
x>0 y>0 z<0
x2/d grand
d grand
*
implicit double precision (a-h)
implicit double precision (o-z)
real coo(3),wg(6)
dimension wgd(6)
parameter (nx=100)
parameter (nmx=2*nx+3)
parameter (eps=1.e-6)
c dimension bjn(0:nmx),bjk (0:nmx)
dimension bjn(0:nmx+nmx),bjk (0:nmx)
dimension a (0:nx) ,bn (0:nx) ,pjn (0:nx), amlt(0:nx)
dimension ax (0:nx) ,bxn (0:nx) ,pjxn(0:nx)
dimension ad (0:nx) ,bdn (0:nx)
dimension bbn (0:nx)
real bek0,bek1,bej0
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
qezs2= 4*exp(z*0.5)
x2sd= x*x*usd
evaluation du nombre de termes
------------------------------
if (x2sd.lt.0.7) then
na=10
elseif (x2sd.lt.3.5) then
na=15
elseif (x2sd.lt.8) then
na=20
elseif (x2sd.lt.14) then
na=25
elseif (x2sd.lt.20) then
na=30
elseif (x2sd.lt.26) then
na=35
else
na=int(x2sd*0.8+25)
endif
evaluation des ji+1/ji(x) [0:nm]
----------------------------------
if (x.lt.5) then
nm=int(4.*X+21.)
elseif (X.lt.10) then
nm=int(2.4*X+28.)
elseif (X.lt.20.) then
nm=int(2.*X+32.)
else
nm=int(1.3*X+52.)
endif
nm=max(nm,na+20)
NM2 = NM+NM
bjn(NM2) =0.
c do 50 i=nm,1,-1
c
c Le nombre de boucles etait faux ; on calculait les termes bjn(i)
c pour i variant de (nm-1) a 0. Or la boucle 60 calcule pjn(i) et pjxn(i)
c pour i variant de 1 a nm en utilisant bjn(2*i). Dans la boucle 220 il y
c a une division par pjxn(i), ce qui, si on allait assez loin dans les indices,
c provoquait une division par 0. (travail de debug en utilisant sigvec et
c un common pour reperer la source du mal...).
c
c ------> La boucle 50 doit donc aller de i=2*nm+1 a 1, et le vecteur bjn
c doit etre declarer deux fois plus grand.
c
do 50 i=NM2+1,1,-1
50 bjn(i-1)=1./( (i+i)/x - bjn(i) )
do 60 i=1,nm
i2=i+i
i2m1=i2-1
pjn(i) =bjn(i2) *bjn(i2m1)
60 pjxn(i) =bjn(i2m1)*bjn(i2-2)
bk0 =dble(bek0(real(ds2)))*exp(-ds2)
bk1 =dble(bek1(real(ds2)))*exp(-ds2)
bj0 =dble(bej0(real(x)))
bj1 =bjn(0)*bj0
bj2 =bjn(1)*bj1
pjn (0) =-1
pjxn(0) =bj0/bj2
recurence montante : an
-----------------------
a(0) =-(qezs2+qezs2)* bj1 *bk0
ax(0) = qezs2 *(bj2-bj0)*bk0
ad(0) = qezs2 * bj1 *bk1
bjk(0)=bk1/bk0
usbjk0=bk0/bk1
bjk(1)= 4*usd+usbjk0
amlt(0) =(pjn(1)-1)*0.5
a (1) =-a (0)*bjk(0)*amlt(0)
ax(1) =-ax(0)*bjk(0)*
* (1+pjxn(1)*(pjxn(2)-2))/(1/pjxn(0)-2+pjxn(1))
ad(1) =-ad(0)*(usbjk0+bjk(1))*0.5*amlt(0)
do 190 i=1,na
ip1=i+1
bjk(ip1)= (ip1+ip1)/ds2+1/bjk(i)
amlt(i) =pjn(i)*(1-pjn(ip1))/(1-pjn(i))
190 continue
naa=na
nax=na
nad=na
do 200 i=1,na
ip1=i+1
a (ip1) =-a (i)*bjk(i)*amlt(i)
if (abs(a(ip1)).lt.eps) then
naa=i
goto 210
endif
200 continue
210 continue
do 220 i=1,na
ip1=i+1
ax(ip1) =-ax(i)*bjk(i)*
* (1+pjxn(ip1)*(pjxn(i+2)-2))/(1/pjxn(i)-2+pjxn(ip1))
if (abs(ax(ip1)).lt.eps) then
nax=i
goto 230
endif
220 continue
230 continue
do 240 i=1,na
ip1=i+1
ad(ip1) =-ad(i)*(1+bjk(ip1)*bjk(i))/(1/bjk(i-1)+bjk(i))*amlt(i)
if (abs(ad(ip1)).lt.eps) then
nad=i
goto 250
endif
240 continue
250 continue
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
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
gwdjk est appelé dans