1 : C KAXK SOURCE CHAT 05/01/13 00:51:57 5004 2 : SUBROUTINE KAXK(A1,A2,OBS,NOBS,NT,NP0,NG0,FF,KIMP,EXTINC,RAD) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C********************************************************************* 6 : C CALCUL DE S1.F12 EN TENANT COMPTE DES OBSTRUCTEURS 7 : C entree 8 : C A1 : COORDONNEES FACE 1 9 : C A2 : COORDONNEES FACE 2 10 : C OBS : COORDONNES DES OBSTRUCTEURS POTENTIELS 11 : C NOBS : NOMBRE D'OBSTRUCTEURS POTENTIELS 12 : C NG0 : NOMBRE DE POINTS DE GAUSS (cas standard) 13 : C NP0 : NOMBRE DE POINTS D'INTEGRATION (elements proches) 14 : C KIMP : parametre 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),OBS(2,NT) 21 : DIMENSION AL(2,20),BL(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 : NM=20 27 : 28 : 29 : C estimation du mode d'integration 30 : NS=2 30 : CALL KAXDIS(A1,A2,NS,KIMP,NG0,NP0,NG,NP) 31 : 32 : IF(KIMP.GE.3) write(6,*) ' kaxk NG NP ',NG,NP 33 : 34 : RI1=A1(1,1) 35 : ZI1=A1(2,1) 36 : RI2=A1(1,2) 37 : ZI2=A1(2,2) 38 : 39 : RJ1=A2(1,1) 40 : ZJ1=A2(2,1) 41 : RJ2=A2(1,2) 42 : ZJ2=A2(2,2) 43 : 44 : DRI=RI2-RI1 45 : DRJ=RJ2-RJ1 46 : DZI=ZI2-ZI1 47 : DZJ=ZJ2-ZJ1 48 : 49 : 50 : C>> MODE D INTEGRATION 51 : IF(NG.EQ.0) THEN 52 : 53 : NA = NP 54 : NB = NP 55 : 56 : C>> INTEGRATION SUR I : A 57 : 58 : FF=0.D0 59 : DA=1./NA 60 : 61 : DO 3 IA=1,NA 62 : 63 : A = DA/2. + DA*(IA-1) 64 : RI=(1.-A)*RI1+A*RI2 65 : ZI=(1.-A)*ZI1+A*ZI2 66 : DA=1./NA 67 : 68 : C>> INTEGRATION SUR J : B 69 : 70 : F=0.D0 71 : DB=1./NB 72 : 73 : DO 30 IB=1,NB 74 : 75 : G=0.D0 76 : B = DB/2. + DB*(IB-1) 77 : RJ=(1.-B)*RJ1+B*RJ2 78 : ZJ=(1.-B)*ZJ1+B*ZJ2 79 : 80 : IF(KIMP.GE.4)WRITE(6,*) ' INTEGRATION IA IB ',IA,IB 81 : 82 : C>> LIMITES DE VISIBILITE PROPRE AUX POINTS I ET J 83 : C ---------------------------------------------- 84 : 84 : CALL KAVOWN(RI,ZI,RJ,ZJ,DRI,DZI,DRJ,DZJ,KVU,NM,NAL,AL,KIMP) 85 : IF(KIMP.GE.4) THEN 86 : WRITE(6,*) ' VISIBILITE PROPRE ',KVU 86 : IF(KVU.NE.0)CALL UTPRIN(AL,2,NAL) 87 : ENDIF 88 : C 89 : C>> TRAITEMENT DES OBSTRUCTEURS 90 : C ---------------------------- 91 : IF(KVU.NE.0.AND.NOBS.NE.0) THEN 91 : CALL KAVOTH(RI,ZI,RJ,ZJ,OBS,NOBS,NT,KVU,NAL,NM,AL,BL,KIMP) 92 : IF(KIMP.GE.4) THEN 93 : WRITE(6,*) ' OBSTRUCTEURS ',KVU 93 : IF(KVU.NE.0)CALL UTPRIN(AL,2,NAL) 94 : ENDIF 95 : ENDIF 96 : C 97 : C>> CALCUL 98 : C ------- 99 : 100 : IF(KVU.NE.0) THEN 100 : CALL KATETA(RI,ZI,RJ,ZJ,DRI,DZI,DRJ,DZJ,NM,NAL,AL,G,KIMP 101 : & ,EXTINC,RAD) 102 : ENDIF 103 : 104 : F= F + 4.*G*DB 105 : IF(KIMP.GE.4)WRITE(6,*) ' IA IB G F ',IA,IB,G,F 106 : 107 : 30 CONTINUE 108 : 109 : FF = FF + F*DA 110 : IF(KIMP.GE.4) WRITE(6,*) ' IA FF ',IA,FF 111 : 112 : 3 CONTINUE 113 : 114 : IF(KIMP.GE.4) WRITE(6,*) ' TOTAL FF ',FF 115 : 116 : ELSE 117 : 118 : C>> POINTS DE GAUSS 119 : 120 : NA = 1 121 : NB = 1 122 : NGA= NG 123 : NGA2=(NGA+1)/2 124 : NGB= NG 125 : NGB2=(NGB+1)/2 125 : CALL MATG(AG) 126 : 127 : IF (AG(1,NGA).LT.1.E-5) THEN 128 : 129 : YA(1)=AG(1,NGA) 130 : HA(1)=AG(2,NGA) 131 : IF(NGA2.GE.2) THEN 132 : DO 100 I=1,NGA2-1 133 : YA(I+1)=AG(2*I+1,NGA) 134 : YA(NGA2+I)=-YA(I+1) 135 : HA(I+1)=AG(2*I+2,NGA) 136 : HA(NGA2+I)=HA(I+1) 137 : 100 CONTINUE 138 : ENDIF 139 : 140 : ELSE 141 : DO 101 I=1,NGA2 142 : YA(I)=AG(2*I-1,NGA) 143 : YA(NGA2+I)=-YA(I) 144 : HA(I)=AG(2*I,NGA) 145 : HA(NGA2+I)=HA(I) 146 : 101 CONTINUE 147 : ENDIF 148 : 149 : IF (AG(1,NGB).LT.1.E-5) THEN 150 : 151 : YB(1)=AG(1,NGB) 152 : HB(1)=AG(2,NGB) 153 : IF(NGB2.GE.2) THEN 154 : DO 200 I=1,NGB2-1 155 : YB(I+1)=AG(2*I+1,NGB) 156 : YB(NGB2+I)=-YA(I+1) 157 : HB(I+1)=AG(2*I+2,NGB) 158 : HB(NGB2+I)=HA(I+1) 159 : 200 CONTINUE 160 : ENDIF 161 : 162 : ELSE 163 : DO 201 I=1,NGB2 164 : YB(I)=AG(2*I-1,NGB) 165 : YB(NGA2+I)=-YB(I) 166 : HB(I)=AG(2*I,NGB) 167 : HB(NGA2+I)=HB(I) 168 : 201 CONTINUE 169 : ENDIF 170 : 171 : C>> INTEGRATION SUR I : A 172 : 173 : FF=0.D0 174 : DA=1./NA 175 : DO 1 IA=1,NA 176 : A = DA/2. + DA*(IA-1) 177 : DA=1./NA 178 : C bornes 179 : AL1=A-DA/2. 180 : AL2=A+DA/2. 181 : 182 : C>> GAUSS SUR I : A 183 : FA=0.D0 184 : DO 11 IGA=1,NGA 185 : C YA varie entre -1 et 1. 186 : ALL= (YA(IGA)+1.)*(AL2-AL1)/2. + AL1 187 : RI=(1.-ALL)*RI1+ALL*RI2 188 : ZI=(1.-ALL)*ZI1+ALL*ZI2 189 : 190 : C>> INTEGRATION SUR J : B 191 : 192 : F=0.D0 193 : DB=1./NB 194 : DO 2 IB=1,NB 195 : B = DB/2. + DB*(IB-1) 196 : C bornes 197 : BL1=B-DB/2. 198 : BL2=B+DB/2. 199 : 200 : C>> GAUSS SUR J : B 201 : FB=0.D0 202 : DO 21 IGB=1,NGB 203 : C YB varie entre -1 et 1. 204 : BLL=(YB(IGB)+1.)*(BL2-BL1)/2. + BL1 205 : RJ=(1.-BLL)*RJ1+BLL*RJ2 206 : ZJ=(1.-BLL)*ZJ1+BLL*ZJ2 207 : 208 : G=0.D0 209 : IF(KIMP.GE.4)WRITE(6,*) ' INTEGRATION IGA IGB ',IGA,IGB 210 : 211 : C>> LIMITES DE VISIBILITE PROPRE AUX POINTS I ET J 212 : C ---------------------------------------------- 213 : 213 : CALL KAVOWN(RI,ZI,RJ,ZJ,DRI,DZI,DRJ,DZJ,KVU,NM,NAL,AL,KIMP) 214 : IF(KIMP.GE.4) THEN 215 : WRITE(6,*) ' VISIBILITE PROPRE ',KVU 215 : IF(KVU.NE.0)CALL UTPRIN(AL,2,NAL) 216 : ENDIF 217 : C 218 : C>> TRAITEMENT DES OBSTRUCTEURS 219 : C ---------------------------- 220 : IF(KVU.NE.0.AND.NOBS.NE.0) THEN 220 : CALL KAVOTH(RI,ZI,RJ,ZJ,OBS,NOBS,NT,KVU,NAL,NM,AL,BL,KIMP) 221 : IF(KIMP.GE.4) THEN 222 : WRITE(6,*) ' OBSTRUCTEURS ',KVU 222 : IF(KVU.NE.0)CALL UTPRIN(AL,2,NAL) 223 : ENDIF 224 : ENDIF 225 : C 226 : C>> CALCUL 227 : C ------- 228 : 229 : IF(KVU.NE.0) THEN 229 : CALL KATETA(RI,ZI,RJ,ZJ,DRI,DZI,DRJ,DZJ,NM,NAL,AL,G,KIMP 230 : & ,EXTINC,RAD) 231 : ENDIF 232 : 233 : FB = FB + 4.*G*HB(IGB)*(BL2-BL1)/2. 234 : 235 : 21 CONTINUE 236 : F= F + FB*DB 237 : 238 : 2 CONTINUE 239 : 240 : FA = FA + F*HA(IGA)*(AL2-AL1)/2. 241 : 242 : 11 CONTINUE 243 : 244 : FF = FF + FA*DA 245 : 246 : 1 CONTINUE 247 : 248 : ENDIF 249 : 250 : RETURN 251 : END 252 : 253 : 254 : 255 :
© Cast3M 2003 - All rights reserved.
Disclaimer