1 : C BANTOU SOURCE PV 09/03/12 21:15:33 6325 2 : SUBROUTINE BANTOU 3 : # (TABCOR,TABVAL,LADIM,ZSIG,VALMOY,DELZET,OMMAX) 4 : C 5 : C----------------------------------------------------------------------- 6 : C Cet opérateur génère un champ aléatoire par la méthode des 7 : C bandes tournantes. 8 : C Renvoi un pointeur actif TABVAL sur la liste des valeurs 9 : C Les coordonnées en entrée ont été adimensionnées par Lambda, 10 : C c-à-d qu'on s'est ramené au cas isotrope avec lambda = 1 11 : C----------------------------------------------------------------------- 12 : C 13 : C--------------------------- 14 : C Paramètres Entrée/Sortie : 15 : C--------------------------- 16 : C 17 : C E/ TABCOR : Pointeur sur la TABLE contenant les coordonnées des 18 : C points où l'on va générer les valeurs (il peut y avoir 19 : C répétition) 20 : C /S TABVAL : Pointeur actif sur la TABLE contenant les valeurs en ces 21 : C points 22 : C E/ LADIM : Indicateur de la dimension relative aux propriétés 23 : C aléatoires du champ à générer 24 : C E/ ZSIG : Ecart-type du champ aléatoire à générer 25 : C E/ DELZET : Incrément d'espace sur la bande (Pas de discrétisation) 26 : C E/ OMMAX : Vecteur d'onde de coupure 27 : C 28 : C-------------------------- 29 : C Variables internes : 30 : C-------------------------- 31 : C 32 : C NLINE : Nombre de bandes 33 : C XLDOM : Dimension maximale du domaine 34 : C DOMEGA : Incrément de vecteur d'onde 35 : C 36 : C--------------------------- 37 : C Ajouts intéressants 38 : C--------------------------- 39 : C 40 : C======================================================================= 41 : C 42 : IMPLICIT INTEGER(I-N) 43 : IMPLICIT REAL*8(A-H,O-Z) 44 : C 45 : -INC CCOPTIO 46 : -INC SMCOORD 47 : -INC SMCHAML 48 : -INC SMELEME 49 : -INC CCREEL 50 : C 51 : C TABLE DES COORDONNEES DES POINTS SUPPORTS 52 : SEGMENT TABCOR 53 : REAL*8 COR(NBL,3) 54 : ENDSEGMENT 55 : C TABLE DES VALEURS EN CES POINTS 56 : SEGMENT TABVAL 57 : REAL*8 VAL(NBL) 58 : ENDSEGMENT 59 : 60 : SEGMENT TH(NLINE) 61 : SEGMENT PHI(NLINE) 62 : SEGMENT ZETMIN(NLINE) 63 : SEGMENT ZETMAX(NLINE) 64 : SEGMENT XLINE(NZET) 65 : C 66 : C EPSILON SERVANT A DIFFERENTS TESTS 67 : C 68 : EPS=1.D-12 69 : C 70 : C INITIALISATION DU CHAMP ALEATOIRE A ZERO 71 : C 72 : SEGACT TABCOR 73 : NBL = COR(/1) 74 : SEGINI TABVAL 75 : DO 20 I=1,NBL 76 : VAL(I) = 0.D0 77 : 20 CONTINUE 78 : C 79 : C =========================================================== 80 : C = Génération des angles décrivant chaque point des lignes = 81 : C =========================================================== 82 : C 83 : C 1D 84 : C ---- 85 : C 86 : IF (LADIM.EQ.1) THEN 87 : C Le champ aléatoire a des propriétés unidirectionnelles 88 : C On génère une seule bande (NLINE = 1) de direction (1 0 0) 89 : C On n'a pas besoin de TH ni de PHI 90 : NLINE=1 91 : C 92 : C 2D 93 : C ---- 94 : C 95 : ELSEIF (LADIM.EQ.2) THEN 96 : C Le champ aléatoire a des propriétés bidimensionnelles 97 : C Générer des lignes de direction répartie régulièrement entre 0 98 : C et 2*PI 99 : C NLINE est figé à une valeur fixe : 100 100 : C 101 : NLINE=100 102 : SEGINI TH 103 : C 104 : C THETA VARIE ENTRE 0 ET 2*PI 105 : C 106 : DO 1 I=1,NLINE 107 : TH(I) = 2.D0*XPI * FLOAT(I-1) / FLOAT(NLINE) 108 : 1 CONTINUE 109 : C 110 : C 3D 111 : C ---- 112 : C 113 : ELSEIF (LADIM.EQ.3) THEN 114 : C Le champ aléatoire a des propriétés tridimensionnelles 115 : C Générer des lignes de direction aléatoire répartie entre 0 et 116 : C 2*PI pour le site, et 0 et PI pour l'azimuth 117 : C NLINE est figé à une valeur fixe : 1000 118 : C 119 : NLINE=1000 120 : SEGINI TH 121 : SEGINI PHI 122 : C 123 : C THETA varie aléatoirement entre 0 et 2*PI, PHI entre 0 et PI 124 : C THETA et COS(PHI) suivent des lois uniformes 125 : C 126 : DO 2 I=1,NLINE 126 : CALL TDRAND(XRAN) 127 : TH(I)=XRAN*XPI*2.D0 127 : CALL TDRAND(XRAN) 128 : PHI(I)=ACOS(1.D0-2.D0*XRAN) 129 : 2 CONTINUE 130 : ENDIF 131 : C 132 : C =========================================================== 133 : C = Détermination de la dimension maximale du domaine XLDOM = 134 : C = et des pas de discrétisation de l'espace dual DOMEGA = 135 : C =========================================================== 136 : C utile uniquement en dimension > 1 137 : C XLDOM est le maximum sur toutes les bandes de l'intervalle des 138 : C valeurs algébriques des positions de la projection de chaque point 139 : C du maillage sur ces droites. 140 : 141 : SEGINI ZETMIN 142 : SEGINI ZETMAX 143 : XLDOM=-1.D0 144 : C 145 : IF (LADIM.EQ.1) THEN 146 : L = 1 147 : ZETMAX(L) = -1 * XGRAND 148 : ZETMIN(L) = XGRAND 149 : DO 9 I=1,NBL 150 : ZET = COR(I,1) 151 : ZETMIN(L) = MIN (ZETMIN(L), ZET) 152 : ZETMAX(L) = MAX (ZETMAX(L), ZET) 153 : 9 CONTINUE 154 : XLDOM = MAX (XLDOM, ZETMAX(L)-ZETMIN(L)) 155 : C 156 : ELSEIF (LADIM.EQ.2) THEN 157 : C 158 : DO 4 L=1,NLINE 159 : ZETMAX(L) = -1 * XGRAND 160 : ZETMIN(L) = XGRAND 161 : THETA = TH(L) 162 : XLI = COS(THETA) 163 : YLI = SIN(THETA) 164 : DO 5 I=1,NBL 165 : X = COR(I,1) 166 : Y = COR(I,2) 167 : ZET = X*XLI + Y*YLI 168 : ZETMIN(L) = MIN (ZETMIN(L), ZET) 169 : ZETMAX(L) = MAX (ZETMAX(L), ZET) 170 : 5 CONTINUE 171 : XLDOM = MAX (XLDOM, ZETMAX(L)-ZETMIN(L)) 172 : 4 CONTINUE 173 : C 174 : ELSE 175 : C 176 : DO 6 L=1,NLINE 177 : ZETMAX(L) = -1 * XGRAND 178 : ZETMIN(L) = XGRAND 179 : THETA = TH(L) 180 : FI = PHI(L) 181 : XLI = COS(THETA) * SIN(FI) 182 : YLI = SIN(THETA) * SIN(FI) 183 : ZLI = COS(FI) 184 : DO 7 I=1,NBL 185 : X = COR(I,1) 186 : Y = COR(I,2) 187 : Z = COR(I,3) 188 : ZET = X*XLI + Y*YLI + Z*ZLI 189 : ZETMIN(L) = MIN (ZETMIN(L), ZET) 190 : ZETMAX(L) = MAX (ZETMAX(L), ZET) 191 : 7 CONTINUE 192 : XLDOM = MAX (XLDOM, ZETMAX(L)-ZETMIN(L)) 193 : 6 CONTINUE 194 : C 195 : ENDIF 196 : C 197 : C DOMEGA doit être inférieur à pi / L 198 : DOMEGA = .9 * (XPI / XLDOM) 199 : C 200 : C =========================================================== 201 : C = Génération du champ aléatoire = 202 : C =========================================================== 203 : C Sommation sur toutes les bandes de la valeur trouvée sur la bande 204 : C en cours au point projeté du point considéré sur cette bande, 205 : C et ce, pour chaque point champ à créer. 206 : C Les valeurs en deux points confondus sont égales. 207 : C 208 : C 1D 209 : C ---- 210 : C 211 : IF (LADIM.EQ.1) THEN 212 : 213 : C Valeurs aléatoires sur la bande 214 : NZET = INT ((ZETMAX(1) - ZETMIN(1)) / DELZET) + 2 215 : SEGINI XLINE 216 : ZETMIL = ZETMIN(1) 216 : CALL BANT1D(DOMEGA,OMMAX,DELZET,ZETMIL,XLINE) 217 : C 218 : C PROJECT LINE PROCESS TO GENERATE 2-D ARRAY 219 : C 220 : DO 8 I=1,NBL 221 : ZET = COR(I,1) 222 : INDEX = INT ((ZET - ZETMIN(1)) / DELZET + 0.5D0) + 1 223 : VAL(I) = VAL(I) + XLINE(INDEX) 224 : 8 CONTINUE 225 : 226 : SEGSUP XLINE 227 : C 228 : C 2D 229 : C ---- 230 : C 231 : ELSEIF (LADIM.EQ.2) THEN 232 : C 233 : DO 12 L=1,NLINE 234 : 235 : C Valeurs aléatoires sur la bande 236 : NZET = INT ((ZETMAX(L) - ZETMIN(L)) / DELZET) + 2 237 : SEGINI XLINE 238 : C 239 : ZETMIL = ZETMIN(L) 239 : CALL BANT2D(DOMEGA,OMMAX,DELZET,ZETMIL,XLINE) 240 : C 241 : C PROJECT LINE PROCESS TO GENERATE 2-D ARRAY 242 : C 243 : THETA = TH(L) 244 : XLI = COS(THETA) 245 : YLI = SIN(THETA) 246 : DO 13 I=1,NBL 247 : X = COR(I,1) 248 : Y = COR(I,2) 249 : ZET = X*XLI + Y*YLI 250 : INDEX = INT ((ZET - ZETMIN(L)) / DELZET + 0.5D0) + 1 251 : VAL(I) = VAL(I) + XLINE(INDEX) 252 : 13 CONTINUE 253 : C 254 : SEGSUP XLINE 255 : C 256 : 12 CONTINUE 257 : 258 : SEGSUP TH 259 : C 260 : C 3D 261 : C ---- 262 : C 263 : ELSEIF (LADIM.EQ.3) THEN 264 : C 265 : DO 15 L=1,NLINE 266 : 267 : C Valeurs aléatoires sur la bande 268 : NZET = INT ((ZETMAX(L) - ZETMIN(L)) / DELZET) + 2 269 : SEGINI XLINE 270 : ZETMIL = ZETMIN(L) 270 : CALL BANT3D(DOMEGA,OMMAX,DELZET,ZETMIL,XLINE) 271 : C 272 : C PROJECT LINE PROCESS TO GENERATE 3-D ARRAY 273 : C 274 : THETA = TH(L) 275 : FI = PHI(L) 276 : XLI = COS(THETA) * SIN(FI) 277 : YLI = SIN(THETA) * SIN(FI) 278 : ZLI = COS(FI) 279 : DO 16 I=1,NBL 280 : X = COR(I,1) 281 : Y = COR(I,2) 282 : Z = COR(I,3) 283 : ZET = X*XLI + Y*YLI + Z*ZLI 284 : INDEX = INT ((ZET - ZETMIN(L)) / DELZET + 0.5D0) + 1 285 : VAL(I) = VAL(I) + XLINE(INDEX) 286 : 16 CONTINUE 287 : C 288 : SEGSUP XLINE 289 : C 290 : 15 CONTINUE 291 : 292 : SEGSUP TH 293 : SEGSUP PHI 294 : ENDIF 295 : C 296 : C.....DIVIDE BY SQRT(FLOAT(NLINE)) 297 : C et mise à la moyenne 298 : C 299 : XX1 = ZSIG / SQRT(FLOAT(NLINE)) 300 : 301 : DO 18 I=1,NBL 302 : VAL(I) = (VAL(I)*XX1) + VALMOY 303 : 18 CONTINUE 304 : C 305 : SEGSUP ZETMIN 306 : SEGSUP ZETMAX 307 : C 308 : RETURN 309 : END 310 : 311 : 312 : 313 : 314 : 315 :
© Cast3M 2003 - All rights reserved.
Disclaimer