[updated 7.May.2002]
SUBROUTINE GHBAND (TYNOYO,K,B,CORM,CORP,E,D1EM,D1EP,D2EMP)
Auteur : O.DeBayser (Septembre 1995)
Derniere modification : D.Martin (7 Mai 2002)
Version 1.0.1
Solution Elementaire au point 'CORM' de l'Equation de HELMHOLTZ
dans une bande de R2 au facteur (i/4) pres:
E(M,P) = Somme d'indice "j" J0(K*Rj)+iY0(K*Rj),
-infini <ou= j <ou= +infini.
-J0 (resp. Y0) est la fonction de Bessel d'ordre 0 de 1e (resp. 2e) espece
-Rj = ( (Xm-Xp)**2 + (Ym-2*j*B-(-1)**j*Yp)**2 )**1/2
-Bande d'axe Ox et de largeur: [-B,+B] en y.
-- Arguments d'entree --
TYNOYO type du noyau de Green a calculer
K constante K de J0(K*Rj)+iY0(K*Rj)
B demi largeur du domaine "bande"
CORM coordonnees du point source
CORP coordonnees du point image
-- Arguments de sortie --
E E(M,P)
D1EM derivees de E par rapport aux coordonnees de M
D1EP derivees de E par rapport aux coordonnees de P
D2EMP derivees secondes croisees dans l'ordre suivant
d2E/dxmdxp, d2E/dxmdyp, d2E/dymdxp, d2E/dymdyp
-- Remarques: La multiplication par i/4 n'est pas effectuee dans ce S-P
On utilise la relation dJ0/dz = -J1 et dY0/dz = -Y1 pour la derivee
On utilise l'equa diff de Bessel : z**2 d2y/dz2 + z dy/dz + z**2 y = 0
pour le calcul de la derivee seconde
CHARACTER*(*) TYNOYO
REAL K,B,CORP(*),CORM(*)
COMPLEX EP,EM,E,D1EP(*),D1EM(*),D2EMP(*)
INTEGER NMAX
REAL MSUN,RJ,X,XJ,AJ0,AJ1,AY0,AY1
COMPLEX EPRIME,ESECON
PARAMETER ( NMAX = 100 )
DATA EPSIL /1.E-10/
E=0.
D1EP(1) = 0.
D1EP(2) = 0.
D1EM(1) = 0.
D1EM(2) = 0.
D2EMP(1)= 0.
D2EMP(2)= 0.
D2EMP(3)= 0.
D2EMP(4)= 0.
MSUN=1.
Pour J=0
X =CORP(1)-CORM(1)
YJ=CORP(2)-CORM(2)
RJ=X*X+YJ*YJ
IF (RJ.LT.1.E-06) GO TO 999
RJ=SQRT(RJ)
UNSUR=1./RJ
Calcul des fonctions de Bessel de 1ere et 2nde espece d'ordre 0 et 1
CALL BEJY01 (K*RJ,AJ0,AJ1,AY0,AY1,.TRUE.)
E=E+CMPLX (AJ0,AY0)
EPRIM = (K/RJ)*(J0prim(K*R) +i YOprim(K*R))
= -(K/RJ)*(J1(K*R) +i Y1(K*R) )
EPRIME=-K*UNSUR*CMPLX(AJ1,AY1)
D1EP(1) = D1EP(1) + X *EPRIME
D1EP(2) = D1EP(2) + YJ*EPRIME
D1EM(1) = D1EM(1) - X *EPRIME
D1EM(2) = D1EM(2) - MSUN*YJ*EPRIME
IF (TYNOYO(:3).EQ.'FOU'.OR.TYNOYO(:3).EQ.'NEU') THEN
ESECON=-K**2 * (J0seconde(K*R) +i YOseconde(K*R))
= K**2 * (J1prim(K*R) +i Y1prim(K*R) )
= K**2(J0(K*R) +i Y0(K*R)) + (K/R)*(J0prim(K*R) +i YOprim(K*R))
= K**2(J0(K*R) +i Y0(K*R)) - (K/R)*(J1(K*R) +i Y1(K*R) )
X =UNSUR*X
YJ=UNSUR*YJ
ESECON=K*K*E +EPRIME+EPRIME
D2EMP(1) =D2EMP(1) + X *X *ESECON- EPRIME
D2EMP(2) =D2EMP(2) + MSUN*YJ*X *ESECON
D2EMP(3) =D2EMP(3) + X *YJ*ESECON
D2EMP(4) =D2EMP(4) + MSUN*YJ*YJ*ESECON-MSUN*EPRIME
ENDIF
DO 50 J=1,NMAX
MSUN=-MSUN
cas indice "j" positif
X =CORP(1)-CORM(1)
YJ=CORP(2)-2*FLOAT(J)*B-MSUN*CORM(2)
RJ=X*X+YJ*YJ
IF(RJ.LT.1.E-06) GO TO 999
RJ=SQRT (RJ)
UNSUR=1./RJ
Calcul des fonctions de Bessel de 1ere et 2nde espece d'ordre 0 et 1
CALL BEJY01 (K*RJ,AJ0,AJ1,AY0,AY1,.TRUE.)
EP=CMPLX(AJ0,AY0)
E=E+EP
EPRIM = (K/RJ)*(J0prim(K*R) +i YOprim(K*R))
= -(K/RJ)*(J1(K*R) +i Y1(K*R) )
EPRIME=-K*UNSUR*CMPLX(AJ1,AY1)
D1EP(1) = D1EP(1) + X *EPRIME
D1EP(2) = D1EP(2) + YJ*EPRIME
D1EM(1) = D1EM(1) - X *EPRIME
D1EM(2) = D1EM(2) - MSUN*YJ*EPRIME
IF (TYNOYO(:3).EQ.'FOU'.OR.TYNOYO(:3).EQ.'NEU') THEN
ESECON=-K**2 * (J0seconde(K*R) +i YOseconde(K*R))
= K**2 * (J1prim(K*R) +i Y1prim(K*R) )
= K**2(J0(K*R) +i Y0(K*R))+ (K/R)*(J0prim(K*R)+i YOprim(K*R))
= K**2(J0(K*R) +i Y0(K*R))- (K/R)*(J1(K*R) +i Y1(K*R))
X =UNSUR*X
YJ=UNSUR*YJ
ESECON=K*K*EP+EPRIME+EPRIME
D2EMP(1) =D2EMP(1) + X *X *ESECON- EPRIME
D2EMP(2) =D2EMP(2) + MSUN*YJ*X *ESECON
D2EMP(3) =D2EMP(3) + X *YJ*ESECON
D2EMP(4) =D2EMP(4) + MSUN*YJ*YJ*ESECON-MSUN*EPRIME
ENDIF
cas indice "j" negatif
X =CORP(1)-CORM(1)
YJ=CORP(2)+2*REAL(J)*B-MSUN*CORM(2)
RJ=X*X+YJ*YJ
IF(RJ.LT.1.E-06) GO TO 999
RJ=SQRT (RJ)
UNSUR=1./RJ
Calcul des fonctions de Bessel de 1ere et 2nde espece d'ordre 0 et 1
CALL BEJY01 (K*RJ,AJ0,AJ1,AY0,AY1,.TRUE.)
EM=CMPLX (AJ0,AY0)
E=E+EM
IF (ABS (EM+EP).LE.EPSIL) GOTO 60
EPRIM = (K/RJ)*(J0prim(K*R) +i YOprim(K*R))
= -(K/RJ)*(J1(K*R) +i Y1(K*R) )
EPRIME=-K*UNSUR*CMPLX(AJ1,AY1)
D1EP(1) = D1EP(1) + X *EPRIME
D1EP(2) = D1EP(2) + YJ*EPRIME
D1EM(1) = D1EM(1) - X *EPRIME
D1EM(2) = D1EM(2) - MSUN*YJ*EPRIME
IF (TYNOYO(:3).EQ.'FOU'.OR.TYNOYO(:3).EQ.'NEU') THEN
ESECON=-K**2 * (J0seconde(K*R) +i YOseconde(K*R))
= K**2 * (J1prim(K*R) +i Y1prim(K*R) )
= K**2(J0(K*R) +i Y0(K*R))+ (K/R)*(J0prim(K*R)+i YOprim(K*R))
= K**2(J0(K*R) +i Y0(K*R))- (K/R)*(J1(K*R) +i Y1(K*R))
X =UNSUR*X
YJ=UNSUR*YJ
ESECON=K*K*EM+EPRIME+EPRIME
D2EMP(1) =D2EMP(1) + X *X *ESECON- EPRIME
D2EMP(2) =D2EMP(2) + MSUN*YJ*X *ESECON
D2EMP(3) =D2EMP(3) + X *YJ*ESECON
D2EMP(4) =D2EMP(4) + MSUN*YJ*YJ*ESECON-MSUN*EPRIME
ENDIF
50 CONTINUE
CALL BAISE ('*Ghband* ---> NON CONVERGENCE de la serie !!!')
STOP
60 CONTINUE
RETURN
999 WRITE (*,1000) (CORP(I),I=1,2),(CORM(I),I=1,2)
CALL BAISE (' Erreur dans le maillage ?')
1000 FORMAT('*GhBand* On essaie de calculer la fonction de Green pour'
&,' une distance < 10**-6'/T11,'XP=',E12.4,' YP=',E12.4,
&' XM=',E12.4,' YM=',E12.4/)
END !Ghband
ghband est appelé dans