1 : C INTFRE SOURCE CHAT 05/01/13 00:41:08 5004 2 : SUBROUTINE INTFRE(CC,SS,XF) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8(A-H,O-Z) 5 : C 6 : C ===================================================================== 7 : C SOUS-PROGRAMME FORTRAN POUR CALCULER LES INTEGRALES DE FRESNEL : 8 : C 9 : C / XF 10 : C CC=1/SQRT(2PI)/ COS(T)/SQRT(T) DT 11 : C / 0 12 : C 13 : C / XF 14 : C SS=1/SQRT(2PI)/ SIN(T)/SQRT(T) DT 15 : C / 0 16 : C 17 : C PROGRAMMEUR : GUILBAUD 18 : C CREATION : 13/01/87 19 : C ===================================================================== 20 : C 21 : -INC CCREEL 22 : DEUPI=2.D0*XPI 23 : PIS2=XPI/2.D0 24 : DX2=XPI/2000.D0 25 : XG=20.D0 26 : XM=10.D0 27 : XP=1.D0 28 : IMAX=10 29 : IF(XF.LE.XM) THEN 30 : XFI=XF 31 : IF(XF.GT.XP) XFI=XP 32 : C 33 : C DEVELOPPEMENT LIMITE POUR XF < XP 34 : C 35 : S2=SQRT(XFI/PIS2) 36 : S0=+XFI/3.D0 37 : SS=S0 38 : C0=+1.D0 39 : CC=C0 40 : XF2=XFI*XFI 41 : DO 10 I=1,IMAX 42 : RS=(XF2*DBLE(4*I-1))/DBLE((4*I+3)*2*I*(2*I+1)) 43 : RC=(XF2*DBLE(4*I-3))/DBLE((4*I+1)*2*I*(2*I-1)) 44 : S0=-S0*RS 45 : C0=-C0*RC 46 : SS=SS+S0 47 : CC=CC+C0 48 : 10 CONTINUE 49 : SS=SS*S2 50 : CC=CC*S2 51 : IF(XF.GT.XP) THEN 52 : C 53 : C INTEGRATION PAR LA METHODE DES TRAPEZES POUR XP < XF < XM 54 : C 55 : KMAX=IDINT((XF-XP)/DX2) 56 : IF(KMAX.EQ.0) KMAX=1 57 : DX=(XF-XP)/DBLE(KMAX+1) 58 : SSX=SIN(XP)*0.5D0/SQRT(XP) 59 : SCX=COS(XP)*0.5D0/SQRT(XP) 60 : XX=XP 61 : DO 20 K=1,KMAX 62 : XX=XX+DX 63 : SSX=SSX+SIN(XX)/SQRT(XX) 64 : SCX=SCX+COS(XX)/SQRT(XX) 65 : 20 CONTINUE 66 : SS=(SSX+0.5D0*SIN(XF)/SQRT(XF))*DX/SQRT(DEUPI)+SS 67 : CC=(SCX+0.5D0*COS(XF)/SQRT(XF))*DX/SQRT(DEUPI)+CC 68 : ENDIF 69 : ELSE 70 : C 71 : C DEVELOPPEMENT ASYMPTOTIQUE POUR XF > XG 72 : C 73 : XFI=XF 74 : IF(XF.LT.XG) XFI=XG 75 : PXX=2.D0*XFI 76 : PXX1=1.D0/PXX 77 : PXX2=3.D0*PXX1/PXX 78 : PXX3=5.D0*PXX2/PXX 79 : PXX4=7.D0*PXX3/PXX 80 : PXX5=9.D0*PXX4/PXX 81 : PXX6=11.D0*PXX5/PXX 82 : PXX7=13.D0*PXX6/PXX 83 : PXX8=15.D0*PXX7/PXX 84 : PXX9=17.D0*PXX8/PXX 85 : PXX10=19.D0*PXX9/PXX 86 : PXX11=21.D0*PXX10/PXX 87 : COEF1=-1.D0+PXX2-PXX4+PXX6-PXX8+PXX10 88 : COEF2=-PXX1+PXX3-PXX5+PXX7-PXX9+PXX11 89 : RPX=SQRT(DEUPI*XFI) 90 : DC=COS(XFI)/RPX 91 : DS=SIN(XFI)/RPX 92 : SS=0.5D0+DC*COEF1+DS*COEF2 93 : CC=0.5D0+DC*COEF2-DS*COEF1 94 : IF(XF.LT.XG) THEN 95 : C 96 : C INTEGRATION PAR LA METHODE DES TRAPEZES POUR XM < XF < XG 97 : C 98 : KMAX=IDINT((XG-XF)/DX2) 99 : IF(KMAX.EQ.0) KMAX=1 100 : DX=(XG-XF)/DBLE(KMAX+1) 101 : SSX=SIN(XF)*0.5D0/SQRT(XF) 102 : SCX=COS(XF)*0.5D0/SQRT(XF) 103 : XX=XF 104 : DO 30 K=1,KMAX 105 : XX=XX+DX 106 : SSX=SSX+SIN(XX)/SQRT(XX) 107 : SCX=SCX+COS(XX)/SQRT(XX) 108 : 30 CONTINUE 109 : SS=SS-(SSX+0.5D0*SIN(XG)/SQRT(XG))*DX/SQRT(DEUPI) 110 : CC=CC-(SCX+0.5D0*COS(XG)/SQRT(XG))*DX/SQRT(DEUPI) 111 : ENDIF 112 : ENDIF 113 : RETURN 114 : END 115 : 116 : 117 :
© Cast3M 2003 - All rights reserved.
Disclaimer