Sources Esope | Cast3M

intfre
C INTFRE    SOURCE    CB215821  17/11/30    21:16:30     9639           
      SUBROUTINE INTFRE(CC,SS,XF)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C
C =====================================================================
C   SOUS-PROGRAMME FORTRAN POUR CALCULER LES INTEGRALES DE FRESNEL :
C
C                  / XF
C   CC=1/SQRT(2PI)/   COS(T)/SQRT(T) DT
C                / 0
C
C                  / XF
C   SS=1/SQRT(2PI)/   SIN(T)/SQRT(T) DT
C                / 0
C
C   PROGRAMMEUR : GUILBAUD
C   CREATION    : 13/01/87
C =====================================================================
C
-INC CCREEL
      DEUPI=2.D0*XPI
      PIS2=XPI/2.D0
      DX2=XPI/2000.D0
      XG=20.D0
      XM=10.D0
      XP=1.D0
      IMAX=10
      IF(XF.LE.XM) THEN
        XFI=XF
        IF(XF.GT.XP) XFI=XP
C
C  DEVELOPPEMENT LIMITE POUR XF &lt; XP
C
        S2=SQRT(XFI/PIS2)
        S0=+XFI/3.D0
        SS=S0
        C0=+1.D0
        CC=C0
        XF2=XFI*XFI
        DO 10 I=1,IMAX
        RS=(XF2*DBLE(4*I-1))/DBLE((4*I+3)*2*I*(2*I+1))
        RC=(XF2*DBLE(4*I-3))/DBLE((4*I+1)*2*I*(2*I-1))
        S0=-S0*RS
        C0=-C0*RC
        SS=SS+S0
        CC=CC+C0
   10   CONTINUE
        SS=SS*S2
        CC=CC*S2
        IF(XF.GT.XP) THEN
C
C  INTEGRATION PAR LA METHODE DES TRAPEZES POUR XP &lt; XF &lt; XM
C
          KMAX=INT((XF-XP)/DX2)
          IF(KMAX.EQ.0) KMAX=1
          DX=(XF-XP)/DBLE(KMAX+1)
          SSX=SIN(XP)*0.5D0/SQRT(XP)
          SCX=COS(XP)*0.5D0/SQRT(XP)
          XX=XP
          DO 20 K=1,KMAX
          XX=XX+DX
          SSX=SSX+SIN(XX)/SQRT(XX)
          SCX=SCX+COS(XX)/SQRT(XX)
  20      CONTINUE
          SS=(SSX+0.5D0*SIN(XF)/SQRT(XF))*DX/SQRT(DEUPI)+SS
          CC=(SCX+0.5D0*COS(XF)/SQRT(XF))*DX/SQRT(DEUPI)+CC
        ENDIF
      ELSE
C
C   DEVELOPPEMENT ASYMPTOTIQUE POUR XF > XG
C
        XFI=XF
        IF(XF.LT.XG) XFI=XG
        PXX=2.D0*XFI
        PXX1=1.D0/PXX
        PXX2=3.D0*PXX1/PXX
        PXX3=5.D0*PXX2/PXX
        PXX4=7.D0*PXX3/PXX
        PXX5=9.D0*PXX4/PXX
        PXX6=11.D0*PXX5/PXX
        PXX7=13.D0*PXX6/PXX
        PXX8=15.D0*PXX7/PXX
        PXX9=17.D0*PXX8/PXX
        PXX10=19.D0*PXX9/PXX
        PXX11=21.D0*PXX10/PXX
        COEF1=-1.D0+PXX2-PXX4+PXX6-PXX8+PXX10
        COEF2=-PXX1+PXX3-PXX5+PXX7-PXX9+PXX11
        RPX=SQRT(DEUPI*XFI)
        DC=COS(XFI)/RPX
        DS=SIN(XFI)/RPX
        SS=0.5D0+DC*COEF1+DS*COEF2
        CC=0.5D0+DC*COEF2-DS*COEF1
        IF(XF.LT.XG) THEN
C
C  INTEGRATION PAR LA METHODE DES TRAPEZES POUR XM &lt; XF &lt; XG
C
          KMAX=INT((XG-XF)/DX2)
          IF(KMAX.EQ.0) KMAX=1
          DX=(XG-XF)/DBLE(KMAX+1)
          SSX=SIN(XF)*0.5D0/SQRT(XF)
          SCX=COS(XF)*0.5D0/SQRT(XF)
          XX=XF
          DO 30 K=1,KMAX
          XX=XX+DX
          SSX=SSX+SIN(XX)/SQRT(XX)
          SCX=SCX+COS(XX)/SQRT(XX)
  30      CONTINUE
          SS=SS-(SSX+0.5D0*SIN(XG)/SQRT(XG))*DX/SQRT(DEUPI)
          CC=CC-(SCX+0.5D0*COS(XG)/SQRT(XG))*DX/SQRT(DEUPI)
        ENDIF
      ENDIF
      RETURN
      END