1 : C KAXC SOURCE CHAT 05/01/13 00:51:51 5004 2 : SUBROUTINE KAXC(A1,A2,NP0,NG0,FF,KIMP,EXTINC,RAD) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C Include contenant quelques constantes dont XPI : 6 : -INC CCREEL 7 : C********************************************************************* 8 : C CALCUL DE S1.F12 EN GEOMETRIE CONVEXE 9 : C entree 10 : C A1 : COORDONNEES FACE 1 11 : C A2 : COORDONNEES FACE 2 12 : C NG0 : NOMBRE DE POINTS DE GAUSS (cas standard) 13 : C NP0 : NOMBRE DE POINTS D INTEGRATION (elements proches) 14 : C KIMP : paramete d'impression 15 : C EXTINC: coefficient d'extinction de la cavite si absorbante 16 : C RAD : dimension du pb (le calcul est fait en coor.reduites) 17 : C resultat 18 : C FF : S1.F12 19 : C********************************************************************* 20 : DIMENSION A1(2,2),A2(2,2) 21 : DIMENSION AL(2,20) 22 : DIMENSION AG(11,10),YA(10),HA(10),YB(10),HB(10) 23 : 24 : C LES INTERVALLES D INTEGRATION SONT AL(I,I+1),I=1,NAL 25 : C 26 : 27 : C estimation du mode d'integration 28 : NS=2 28 : CALL KAXDIS(A1,A2,NS,KIMP,NG0,NP0,NG,NP) 29 : 30 : IF(KIMP.GE.3) write(6,*) ' kaxc NG NP ',NG,NP 31 : 32 : 33 : NM=20 34 : NAL=1 35 : AL(1,1) = 0.D0 36 : AL(2,1) = XPI 37 : 38 : RI1=A1(1,1) 39 : ZI1=A1(2,1) 40 : RI2=A1(1,2) 41 : ZI2=A1(2,2) 42 : 43 : RJ1=A2(1,1) 44 : ZJ1=A2(2,1) 45 : RJ2=A2(1,2) 46 : ZJ2=A2(2,2) 47 : 48 : DRI=RI2-RI1 49 : DRJ=RJ2-RJ1 50 : DZI=ZI2-ZI1 51 : DZJ=ZJ2-ZJ1 52 : 53 : 54 : C>> MODE D INTEGRATION 55 : IF(NG.EQ.0) THEN 56 : 57 : NA = NP 58 : NB = NP 59 : 60 : C>> INTEGRATION SUR I : A 61 : 62 : FF=0.D0 63 : DA=1.D0/NA 64 : 65 : DO 3 IA=1,NA 66 : 67 : A = DA/2.D0 + DA*(IA-1) 68 : RI=(1.D0-A)*RI1+A*RI2 69 : ZI=(1.D0-A)*ZI1+A*ZI2 70 : DA=1.D0/NA 71 : 72 : C>> INTEGRATION SUR J : B 73 : 74 : F=0.D0 75 : DB=1.D0/NB 76 : 77 : DO 30 IB=1,NB 78 : 79 : G=0.D0 80 : B = DB/2.D0 + DB*(IB-1) 81 : RJ=(1.D0-B)*RJ1+B*RJ2 82 : ZJ=(1.D0-B)*ZJ1+B*ZJ2 83 : 84 : IF(KIMP.GE.4)WRITE(6,*) ' INTEGRATION IA IB ',IA,IB 85 : C 86 : C>> CALCUL 87 : C ------- 88 : 88 : CALL KATETA(RI,ZI,RJ,ZJ,DRI,DZI,DRJ,DZJ,NM,NAL,AL,G,KIMP 89 : & ,EXTINC,RAD) 90 : 91 : F= F + 4.D0*G*DB 92 : IF(KIMP.GE.4)WRITE(6,*) ' IA IB G F ',IA,IB,G,F 93 : 94 : 30 CONTINUE 95 : 96 : FF = FF + F*DA 97 : IF(KIMP.GE.4) WRITE(6,*) ' IA FF ',IA,FF 98 : 99 : 3 CONTINUE 100 : 101 : IF(KIMP.GE.4) WRITE(6,*) ' TOTAL FF ',FF 102 : 103 : ELSE 104 : 105 : C>> POINTS DE GAUSS 106 : 107 : NA = 1 108 : NB = 1 109 : NGA= NG 110 : NGA2=(NGA+1)/2 111 : NGB= NG 112 : NGB2=(NGB+1)/2 112 : CALL MATG(AG) 113 : 114 : IF (AG(1,NGA).LT.1.D-5) THEN 115 : 116 : YA(1)=AG(1,NGA) 117 : HA(1)=AG(2,NGA) 118 : IF(NGA2.GE.2) THEN 119 : DO 100 I=1,NGA2-1 120 : YA(I+1)=AG(2*I+1,NGA) 121 : YA(NGA2+I)=-YA(I+1) 122 : HA(I+1)=AG(2*I+2,NGA) 123 : HA(NGA2+I)=HA(I+1) 124 : 100 CONTINUE 125 : ENDIF 126 : 127 : ELSE 128 : DO 101 I=1,NGA2 129 : YA(I)=AG(2*I-1,NGA) 130 : YA(NGA2+I)=-YA(I) 131 : HA(I)=AG(2*I,NGA) 132 : HA(NGA2+I)=HA(I) 133 : 101 CONTINUE 134 : ENDIF 135 : 136 : IF (AG(1,NGB).LT.1.D-5) THEN 137 : 138 : YB(1)=AG(1,NGB) 139 : HB(1)=AG(2,NGB) 140 : IF(NGB2.GE.2) THEN 141 : DO 200 I=1,NGB2-1 142 : YB(I+1)=AG(2*I+1,NGB) 143 : YB(NGB2+I)=-YA(I+1) 144 : HB(I+1)=AG(2*I+2,NGB) 145 : HB(NGB2+I)=HA(I+1) 146 : 200 CONTINUE 147 : ENDIF 148 : 149 : ELSE 150 : DO 201 I=1,NGB2 151 : YB(I)=AG(2*I-1,NGB) 152 : YB(NGA2+I)=-YB(I) 153 : HB(I)=AG(2*I,NGB) 154 : HB(NGA2+I)=HB(I) 155 : 201 CONTINUE 156 : ENDIF 157 : 158 : C>> INTEGRATION SUR I : A 159 : 160 : FF=0.D0 161 : DA=1.D0/NA 162 : DO 1 IA=1,NA 163 : A = DA/2.D0 + DA*(IA-1) 164 : DA=1.D0/NA 165 : C bornes 166 : AL1=A-DA/2.D0 167 : AL2=A+DA/2.D0 168 : 169 : C>> GAUSS SUR I : A 170 : FA=0.D0 171 : DO 11 IGA=1,NGA 172 : C YA varie entre -1 et 1.D0 173 : ALL= (YA(IGA)+1.D0)*(AL2-AL1)/2.D0 + AL1 174 : RI=(1.D0-ALL)*RI1+ALL*RI2 175 : ZI=(1.D0-ALL)*ZI1+ALL*ZI2 176 : 177 : C>> INTEGRATION SUR J : B 178 : 179 : F=0.D0 180 : DB=1.D0/NB 181 : DO 2 IB=1,NB 182 : B = DB/2.D0 + DB*(IB-1) 183 : C bornes 184 : BL1=B-DB/2.D0 185 : BL2=B+DB/2.D0 186 : 187 : C>> GAUSS SUR J : B 188 : FB=0.D0 189 : DO 21 IGB=1,NGB 190 : C YB varie entre -1 et 1.D0 191 : BLL=(YB(IGB)+1.D0)*(BL2-BL1)/2.D0 + BL1 192 : RJ=(1.D0-BLL)*RJ1+BLL*RJ2 193 : ZJ=(1.D0-BLL)*ZJ1+BLL*ZJ2 194 : 195 : G=0.D0 196 : C 197 : C>> CALCUL 198 : C ------- 199 : 199 : CALL KATETA(RI,ZI,RJ,ZJ,DRI,DZI,DRJ,DZJ,NM,NAL,AL,G 200 : $ ,KIMP,EXTINC,RAD) 201 : 202 : FB = FB + 4.D0*G*HB(IGB)*(BL2-BL1)/2.D0 203 : 204 : 21 CONTINUE 205 : F= F + FB*DB 206 : 207 : 2 CONTINUE 208 : 209 : FA = FA + F*HA(IGA)*(AL2-AL1)/2.D0 210 : 211 : 11 CONTINUE 212 : 213 : FF = FF + FA*DA 214 : 215 : 1 CONTINUE 216 : 217 : ENDIF 218 : 219 : RETURN 220 : END 221 : 222 : 223 : 224 : 225 : 226 : 227 : 228 :
© Cast3M 2003 - All rights reserved.
Disclaimer