Sources Esope | Cast3M

kaxc
C KAXC      SOURCE    CHAT      05/01/13    00:51:51     5004
      SUBROUTINE KAXC(A1,A2,NP0,NG0,FF,KIMP,EXTINC,RAD)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C    Include contenant quelques constantes dont XPI :
-INC CCREEL
C*********************************************************************
C     CALCUL DE  S1.F12   EN GEOMETRIE CONVEXE
C    entree
C     A1    :  COORDONNEES FACE 1
C     A2    :  COORDONNEES FACE 2
C     NG0   :  NOMBRE DE POINTS DE GAUSS (cas standard)
C     NP0   :  NOMBRE DE POINTS D INTEGRATION (elements proches)
C     KIMP  :  paramete d'impression
C     EXTINC:  coefficient d'extinction de la cavite si absorbante
C     RAD   :  dimension du pb (le calcul est fait en coor.reduites)
C    resultat
C     FF    : S1.F12
C*********************************************************************
      DIMENSION A1(2,2),A2(2,2)
      DIMENSION AL(2,20)
      DIMENSION AG(11,10),YA(10),HA(10),YB(10),HB(10)
 
C    LES INTERVALLES D INTEGRATION SONT AL(I,I+1),I=1,NAL
C
 
C    estimation du mode d'integration
      NS=2
      CALL KAXDIS(A1,A2,NS,KIMP,NG0,NP0,NG,NP)
 
      IF(KIMP.GE.3) write(6,*) ' kaxc NG NP ',NG,NP
 
 
      NM=20
      NAL=1
      AL(1,1) = 0.D0
      AL(2,1) = XPI
 
      RI1=A1(1,1)
      ZI1=A1(2,1)
      RI2=A1(1,2)
      ZI2=A1(2,2)
 
      RJ1=A2(1,1)
      ZJ1=A2(2,1)
      RJ2=A2(1,2)
      ZJ2=A2(2,2)
 
      DRI=RI2-RI1
      DRJ=RJ2-RJ1
      DZI=ZI2-ZI1
      DZJ=ZJ2-ZJ1
 
 
C>> MODE D INTEGRATION
      IF(NG.EQ.0) THEN
 
         NA = NP
         NB = NP
 
C>>   INTEGRATION SUR I : A
 
         FF=0.D0
         DA=1.D0/NA
 
         DO 3  IA=1,NA
 
            A = DA/2.D0 + DA*(IA-1)
            RI=(1.D0-A)*RI1+A*RI2
            ZI=(1.D0-A)*ZI1+A*ZI2
            DA=1.D0/NA
 
C>>   INTEGRATION SUR J : B
 
            F=0.D0
            DB=1.D0/NB
 
            DO 30 IB=1,NB
 
               G=0.D0
               B = DB/2.D0 + DB*(IB-1)
               RJ=(1.D0-B)*RJ1+B*RJ2
               ZJ=(1.D0-B)*ZJ1+B*ZJ2
 
               IF(KIMP.GE.4)WRITE(6,*) ' INTEGRATION IA IB ',IA,IB
C
C>>   CALCUL
C    -------
 
               CALL KATETA(RI,ZI,RJ,ZJ,DRI,DZI,DRJ,DZJ,NM,NAL,AL,G,KIMP
     &              ,EXTINC,RAD)
 
               F= F + 4.D0*G*DB
               IF(KIMP.GE.4)WRITE(6,*) ' IA IB G F ',IA,IB,G,F
 
 30         CONTINUE
 
            FF = FF + F*DA
            IF(KIMP.GE.4) WRITE(6,*) ' IA  FF ',IA,FF
 
 3       CONTINUE
 
         IF(KIMP.GE.4) WRITE(6,*) ' TOTAL FF ',FF
 
      ELSE
 
C>> POINTS DE GAUSS
 
         NA = 1
         NB = 1
         NGA= NG
         NGA2=(NGA+1)/2
         NGB= NG
         NGB2=(NGB+1)/2
         CALL MATG(AG)
 
         IF (AG(1,NGA).LT.1.D-5) THEN
 
            YA(1)=AG(1,NGA)
            HA(1)=AG(2,NGA)
            IF(NGA2.GE.2) THEN
               DO 100 I=1,NGA2-1
                  YA(I+1)=AG(2*I+1,NGA)
                  YA(NGA2+I)=-YA(I+1)
                  HA(I+1)=AG(2*I+2,NGA)
                  HA(NGA2+I)=HA(I+1)
 100           CONTINUE
            ENDIF
 
         ELSE
            DO 101 I=1,NGA2
               YA(I)=AG(2*I-1,NGA)
               YA(NGA2+I)=-YA(I)
               HA(I)=AG(2*I,NGA)
               HA(NGA2+I)=HA(I)
 101        CONTINUE
         ENDIF
 
         IF (AG(1,NGB).LT.1.D-5) THEN
 
            YB(1)=AG(1,NGB)
            HB(1)=AG(2,NGB)
            IF(NGB2.GE.2) THEN
               DO 200 I=1,NGB2-1
                  YB(I+1)=AG(2*I+1,NGB)
                  YB(NGB2+I)=-YA(I+1)
                  HB(I+1)=AG(2*I+2,NGB)
                  HB(NGB2+I)=HA(I+1)
 200           CONTINUE
            ENDIF
 
         ELSE
            DO 201 I=1,NGB2
               YB(I)=AG(2*I-1,NGB)
               YB(NGA2+I)=-YB(I)
               HB(I)=AG(2*I,NGB)
               HB(NGA2+I)=HB(I)
 201        CONTINUE
         ENDIF
 
C>>   INTEGRATION SUR I : A
 
         FF=0.D0
         DA=1.D0/NA
         DO 1  IA=1,NA
            A = DA/2.D0 + DA*(IA-1)
            DA=1.D0/NA
C bornes
            AL1=A-DA/2.D0
            AL2=A+DA/2.D0
 
C>>   GAUSS  SUR I : A
            FA=0.D0
            DO 11 IGA=1,NGA
C YA varie entre -1 et 1.D0
               ALL= (YA(IGA)+1.D0)*(AL2-AL1)/2.D0 + AL1
               RI=(1.D0-ALL)*RI1+ALL*RI2
               ZI=(1.D0-ALL)*ZI1+ALL*ZI2
 
C>>   INTEGRATION SUR J : B
 
               F=0.D0
               DB=1.D0/NB
               DO 2  IB=1,NB
                  B = DB/2.D0 + DB*(IB-1)
C bornes
                  BL1=B-DB/2.D0
                  BL2=B+DB/2.D0
 
C>>   GAUSS  SUR J : B
                  FB=0.D0
                  DO 21 IGB=1,NGB
C YB varie entre -1 et 1.D0
                     BLL=(YB(IGB)+1.D0)*(BL2-BL1)/2.D0 + BL1
                     RJ=(1.D0-BLL)*RJ1+BLL*RJ2
                     ZJ=(1.D0-BLL)*ZJ1+BLL*ZJ2
 
                     G=0.D0
C
C>>   CALCUL
C    -------
 
                     CALL KATETA(RI,ZI,RJ,ZJ,DRI,DZI,DRJ,DZJ,NM,NAL,AL,G
     $                    ,KIMP,EXTINC,RAD)
 
                     FB = FB + 4.D0*G*HB(IGB)*(BL2-BL1)/2.D0
 
 21               CONTINUE
                  F= F + FB*DB
 
 2             CONTINUE
 
               FA = FA + F*HA(IGA)*(AL2-AL1)/2.D0
 
 11         CONTINUE
 
            FF = FF + FA*DA
 
 1       CONTINUE
 
      ENDIF
 
      RETURN
      END