1 : C SIMUL1 SOURCE BP208322 12/01/10 21:15:08 7244 2 : SUBROUTINE SIMUL1 (FREQ,NBMODE,IPRIGI,IPMASS, IPSOLU, LIMAGE) 3 : 4 : ************************************************************************ 5 : * 6 : * S I M U L 1 7 : * ----------- 8 : * 9 : * FONCTION: 10 : * --------- 11 : * 12 : * APPELE PAR LE SOUS-PROGRAMME "SIMULT". 13 : * DETERMINE UNE SERIE DE MODES PROPRES DONT LES FREQUENCES SONT 14 : * VOISINES D'UNE VALEUR DONNEE. 15 : * 16 : * MODE D'APPEL: 17 : * ------------- 18 : * 18 : * CALL SIMUL1 (FREQ,NBMODE,IPRIGI,IPMASS, IPSOLU) 19 : * 20 : * PARAMETRES: (E)=ENTREE (S)=SORTIE 21 : * ----------- 22 : * 23 : * FREQ REEL DP (E) FREQUENCE AUTOUR DE LAQUELLE ON CHERCHE 24 : * DES FREQUENCES PROPRES. 25 : * NBMODE ENTIER (E) NOMBRE DE MODES PROPRES DEMANDES. 26 : * IPRIGI ENTIER (E) POINTEUR SUR L'OBJET 'RIGIDITE' DE 27 : * SOUS-TYPE "RIGIDITE". 28 : * IPMASS ENTIER (E) POINTEUR SUR L'OBJET 'RIGIDITE' DE 29 : * SOUS-TYPE "MASSE". 30 : * IPSOLU ENTIER (S) POINTEUR SUR L'OBJET 'SOLUTION' CONTENANT 31 : * LES MODES PROPRES CALCULES. 32 : * LIMAGE (E) LOGICAL =vrai si on souhaite garder les eventuelles 33 : * valeurs propres negatives 34 : * (matrice K non definie positive par ex.) 35 : * 36 : * MODE DE FONCTIONNEMENT: 37 : * ----------------------- 38 : * 39 : * . SELON LA MEMOIRE, ON SE DONNE UN NOMBRE MAX DE VECTEURS X 40 : * . PAR CYCLE, ON REALISE LES OPERATIONS SUIVANTES : 41 : * 1. (DECALE) DECALAGE DE LA MATRICE DE RIGIDITE K'= K -(2*pi*freq)^2*M 42 : * 2. (LANCZO) CALCUL DES VECTEURS K'-ORTHOGONAUX DE LANCZOS = X 43 : * (eventuellement M-orthogonaux a qq Z deja convergés) 44 : * 3. (SIMU21) CALCUL DES VALEURS et VECTEURS PROPRES DU PROBLEME REDUIT 45 : * 4. (SIMUL3) DEDUCTION DES FREQUENCES PROPRES DU PROBLEME PHYSIQUE, 46 : * RECOMBINAISON DES VECTEURS PROPRES Z=X*PHI, 47 : * CALCUL DU RESIDU, 48 : * ORTHOGONALISATION % MODES DEJA CONVERGES 49 : * 5. (SIMUL6) CREATION DE L'OBJET "SOLUTION" DE CE CYCLE 50 : * 6. (SIMUL7) CALCUL DES MASSE ET DEP GENERALISES ET NUMERO DE MODES 51 : * 7. (STRATE) ADJONCTION DES MODES A CEUX DES CYCLES PRECEDENTS, 52 : * CALCUL D'UN SHIFT A APPLIQUER AU PROCHAIN CYCLE, 53 : * VERIFICATION DE LA COMPLETUDE DU SPECTRE 54 : * 55 : * DESCRIPTION DES VARIABLES : cf. commentaires au cours du programme 56 : * -------------------------- 57 : * 58 : * CREATION, MODIFICATIONS : 59 : * ------------------------ 60 : * 61 : * CREATION : PASCAL MANIGOT 22 AVRIL 1985 62 : * REFONTE : THIERRY CHARRAS, BENOIT PRABEL 2010 63 : * 64 : * LANGAGE : ESOPE (FORTRAN77) 65 : * 66 : ************************************************************************ 67 : * 68 : IMPLICIT INTEGER(I-N) 69 : IMPLICIT REAL*8 (A-H,O-Z) 70 : -INC CCREEL 71 : -INC CCOPTIO 72 : -INC SMLMOTS 73 : -INC SMRIGID 74 : -INC SMTABLE 75 : -INC SMVECTD 76 : -INC SMCHPOI 77 : -INC SMLENTI 78 : c -INC SMLCHPO 79 : c -INC SMSOLUT 80 : -INC SMMATRI 81 : * 82 : * REGOUPEMENT DE RENSEIGNEMENTS POUR "LANCZO" et "SIMUL.": 83 : COMMON/CLANCZ/ IPV1,IPMV1,IPLMOX,IPLMOY 84 : COMMON/CSIMUL/ W2SHIF,XPREC21,mvecri,IPVECX,IPVECZ,IPMZ,IPFREZ 85 : * 86 : REAL*8 W2SHIF,FREQ,DEUXPI,XPREC21 87 : LOGICAL LIMAGE,LIMAG2 88 : LOGICAL lshift 89 : INTEGER OOOVAL 90 : * 91 : PARAMETER (DEUXPI = (2.D0*XPI)) 92 : * 93 : * 94 : * -- INITIALISATIONS -- 95 : * 96 : * du temps cpu passé dans SIMUL1 et les subroutines appelés 96 : call gibtem(xkt) 97 : xkt1 = 0.D0 98 : 99 : * de la precision demandée a simu21 a simul3 100 : XPREC21 = EPSLON(1.0D0) 101 : XPREC21 = XPREC21**0.5 102 : XPREC21 = max(XPREC21,1.D-8) 103 : 104 : * liste des lambdas, des vect p, des permutations du cycle en cours 105 : IPVALP = 0 106 : IPHI = 0 107 : IPERM = 0 108 : * nombre de vect p converges pour le cycle / le shift en cours 109 : nbonve = 0 110 : nbonZ = 0 111 : * solution du cycle en cours / cumul 112 : IPSOLU = 0 113 : idist = 0 114 : * cumul des vect p converges Z, et produit M*Z / |Z^T*M*Z| 115 : IPVECZ = 0 116 : IPMZ = 0 117 : * nombre total de vect p ayant deja converges 118 : nvecZ = 0 119 : ifini = 0 120 : 121 : * decalage spectral initial =shifti et courant=FREQ(->W2SHIF) 122 : lshift = .false. 123 : W2SHIF = (DEUXPI * FREQ) ** 2 124 : IF (LIMAGE) THEN 125 : W2SHIF = SIGN(W2SHIF,FREQ) 126 : ENDIF 127 : shifti = freq 128 : frshif = shifti 128 : CALL DECALE (IPRIGI,IPMASS,W2SHIF, IPKW2M) 129 : IPKWM0= IPKW2M 130 : 131 : * vecteur aleatoire V et M*V 131 : call aleat1(IPKW2M,IPV1) 131 : CALL MUCPRI(IPV1,IPMASS,IPMV1) 131 : CALL CORRSP(IPKW2M,IPV1,IPMV1,IPLMOX,IPLMOY) 131 : CALL INITFL(IPKW2M,IPMASS,IPMV1,IPV1,IFLU) 132 : 133 : *ajout bp 06/01/2012: 134 : * M est symetriques, mais est elle définie positive? 134 : CALL DIAGN1(IPMASS,nvp0M) 135 : * correction pour elements fluides (inconnue PI mise à 0 via INITFL) 136 : c nvp0M = nvp0M - IFLU 137 : * bp 10/01/2012: nvp0M et NEMSM ne semblent pas bien calculés ... 138 : * (resultats dependant machine -> cf. dyna7.dgibi) 139 : * on propose la solution qui suppose que nvp0M si M est LIQUIDE 140 : if(IFLU.gt.0) nvp0M=0 141 : * M est elle singuliere ? 142 : MRIGID=IPMASS 143 : segact,MRIGID 144 : MMATRI=ICHOLE 145 : segdes,MRIGID 146 : segact,MMATRI 147 : NENSM = NENS 148 : segdes,MMATRI 149 : * K est symetriques, mais est elle définie positive? 149 : CALL DIAGN1(IPRIGI,nvp0K) 150 : if(iimpi.ge.5) write(6,*) 'nbre de vp<0 pour M (',IPMASS,') =', 151 : & nvp0M,' et K (',IPRIGI,') =',nvp0K 152 : IF (nvp0M.ne.0.and.nvp0K.ne.0) THEN 153 : write(6,*) 'AU MOINS L UNE DES 2 MATRICES (MASSE OU RIGIDITE)' 154 : write(6,*) 'DOIT ETRE DEFINIE POSITIVE POUR GARANTIR UN BON ' 155 : & ,'FONCTIONNEMENT DE VIBR SIMUL !' 156 : write(6,*) 'L EXECUTION CONTINUE BIEN QUE ' 157 : & ,'LA NUMEROTATION CORRECTE SOIT IMPOSSIBLE...' 158 : c on aurait pu faire une erreur ici car la numerotation est vitale 159 : c dans strate et le numero donné par simul7 est necessairement faux... 159 : c call ERREUR(21) 160 : ENDIF 161 : 162 : *ajout bp 13/05/2011: on laisse simul3 calculer toutes les val p >0 et <0. 163 : * strate fera eventuellement le tri (si .not.limage) 164 : * en fonction de nvp0 (= numero de la dernière val p negative) 165 : LIMAG2=.true. 166 : IF (LIMAGE) THEN 167 : nvp0 = 0 168 : ELSE 169 : if (nvp0M.eq.0) then 170 : nvp0 = nvp0K 171 : elseif(nvp0K.eq.0) then 172 : nvp0 = nvp0M 173 : else 174 : c dans ce cas on met nvp0 a 0 175 : nvp0 = 0 176 : endif 177 : * => il y a donc nvp0 valeurs propres lambda < 0 178 : * on inclut les mode de corps rigides (tq lambda=0) 179 : ENDIF 180 : c if(iimpi.ge.4) write(6,*) '=> nvp0=',nvp0 181 : 182 : * besoin de M-orthogonaliser par rapport aux modes deja convergés ? 183 : Northo = 0 184 : 185 : * TCPU pour decomposer K^shift = LDL^T et creer V aleat 186 : xktale=0.D0 186 : call gibtem(xktale) 187 : xkt1 = xkt1 + xktale 188 : xktal1 = xktale 189 : 190 : * maxmem=taille memoire , inc=nombre d'inconnues du pb 191 : maxmem=oooval(1,1) 192 : mchpoi=IPV1 193 : segact mchpoi 194 : inc=0 195 : do iou=1,ipchp(/1) 196 : msoupo=ipchp(iou) 197 : segact msoupo 198 : mpoval=ipoval 199 : segact mpoval 200 : inc=inc + vpocha(/1)*vpocha(/2) 201 : segdes mpoval,msoupo 202 : enddo 203 : segdes mchpoi 204 : 205 : * incvrai = nombre d'inconnues du pb independantes 206 : * (recopié depuis chole.eso et simplifié..?) 207 : MRIGID=IPKW2M 208 : segact,MRIGID 209 : MMATRI=ICHOLE 210 : segdes,MRIGID 211 : segact,MMATRI 212 : MILIGN=IILIGN 213 : segdes,MMATRI 214 : segact,MILIGN 215 : NBLAG=0 216 : c on espere qu'il faut bien regarder tout le tableau ITTR 217 : c (boucle + compliqué dans chole.eso) 218 : NITTR=ITTR(/1) 219 : do II=1,NITTR 220 : IF(ITTR(II).NE.0) NBLAG=NBLAG+1 221 : enddo 222 : segdes,MILIGN 223 : incvrai = inc - (3*NBLAG/2) 224 : c il faut egalement deduire le nombre d'inconnu PI 225 : incvrai = incvrai - IFLU 226 : c ainsi que le nombre NENSM de valeur propre infinie (M singuliere) 227 : c (rem bp : cela ne semble pas suffire au calcul de incvrai 228 : c NENSM n'est pas toujours exact...? ou mauvais raisonnement...?) 229 : incvrai = incvrai - NENSM 230 : if(iimpi.ge.5) write(6,*)'incvrai=inc-1.5*NBLAG-IFLU-NENSM' 231 : & ,incvrai,inc,NBLAG,IFLU,NENSM 232 : 233 : * Nbmode = nombre de modes recherchés au total 234 : c Nbmode=min(Nbmode,incvrai) 235 : if (incvrai.lt.Nbmode) then 236 : write(6,*) ' simul1 : nombre de modes recherches ', 237 : $ '> nombre d inconnues independantes' 238 : Nbmode = incvrai 239 : write(6,*) ' on reduit le nombre de modes recherches a ',Nbmode 240 : endif 241 : * nvecma = nombre de vecteurs X max. % memoire 242 : Nvecma= maxmem /inc / 30 243 : * nombre de cycle a réaliser 244 : nfois= Nbmode / Nvecma +1 245 : nfois4 = nfois * 1000 246 : cbp: 32*Nbmode car pb lorsque Nbmode=1 et mode multiple! 247 : nfois4 = min(nfois4,32*Nbmode) 248 : Nvecma= min(inc,Nvecma) 249 : * Nmod = nombre de modes recherchés / cycle 250 : Nmod= min(Nbmode,Nvecma/2) 251 : Nbmode0=Nbmode 252 : * Nmodopt = nombre de mode optimal % efficacité tcpu, Nmodma = limite 253 : Nmodopt = Nmod 254 : Nmodma = Nmod 255 : * Nmodcyc = nombre de modes / cycle demandés par l utilisateur 255 : CALL LIRENT(Nmodcyc,0,IRETOU) 256 : if (IRETOU.ne.0) then 257 : Nmod = Nmodcyc 258 : Nmodopt = Nmodcyc 259 : endif 260 : 261 : if (iimpi.ge.4) then 262 : write(6,*) '=============================================' 263 : write(6,*) ' VIBR SIMUL : METHODE DE LANCZOS avec Cycles ' 264 : write(6,*) '=============================================' 265 : write(6,*) ' taille memoire , nb inconnues (vraies)', 266 : $ maxmem,inc,' (',incvrai,')' 267 : write(6,*) ' nbre maxi de vecteurs en simultanée ' , nvecma 268 : write(6,*) ' nbre modes recherchés et /cycle ',Nbmode,Nmod 269 : write(6,*) ' nbre de cycle mini , maxi ' , nfois,nfois4 270 : write(6,*) '=============================================' 271 : write(6,*) ' icycle | Nbre de modes freq ', 272 : $ ' | Nbre de modes | Nbre de modes ' 273 : write(6,*) ' icycle | cherchés shift ', 274 : $ ' | trouvés | total' 275 : endif 276 : if(iimpi.ge.6) write(6,*) 'precision demandee a simu21=',XPREC21 277 : * 278 : * 279 : *-----Cycle de Lanczos------------------------------------------------* 280 : * 281 : do icycle=1,nfois4 282 : 283 : if (iimpi.ge.5) then 284 : write(6,*) '__________________________________________________' 285 : write(6,*) 'Simul1: icycle nmod W2SHIF', icycle, Nmod, W2SHIF 286 : endif 287 : 288 : * 289 : * ------------------------------ 290 : * -- INITIALISATIONS DIVERSES -- 291 : * 292 : * decalage spectral + vecteur initial aleatoire 293 : if (icycle.ne.1) then 294 : cbp10/2010 : ipsolu mis a 0 a tous les cycles (pas seulement si on shifte) 295 : c mais nbonZ = nbre de mode pour ce shift (=plusieurs cycles) 296 : c nbonZ = 0 297 : IPSOLU = 0 298 : * strate peut avoir changé Nbmode et IPKW2M 299 : Nbmode = Nbmode0 300 : if (lshift) then 300 : CALL DECALE(IPRIGI,IPMASS,W2SHIF, IPKW2M) 301 : IPKWM0= IPKW2M 302 : lshift = .false. 303 : nbonZ = 0 304 : c IPSOLU = 0 305 : * l utilisateur 'force' le nombre de modes / cycle 306 : if (IRETOU.ne.0) then 307 : Nmod = Nmodcyc 308 : Nmodopt = Nmodcyc 309 : endif 310 : endif 310 : call aleat1(IPKW2M,IPV1) 310 : CALL MUCPRI (IPV1,IPMASS, IPMV1) 310 : CALL CORRSP(IPKW2M,IPV1,IPMV1,IPLMOX,IPLMOY) 310 : CALL INITFL(IPKW2M,IPMASS,IPMV1,IPV1,IFLU) 311 : endif 312 : IF(IERR.NE.0) RETURN 313 : * 314 : * verrue pour le cas ou l utilisateur force le nombre de modes/cycle 315 : icyc1 = icycle 316 : if(IRETOU.ne.0) icyc1=1 317 : 318 : * 319 : * ce sous programme na plus raison d etre: a supprimer prochainement 319 : * CALL SIMUL0 320 : * 321 : 322 : * ------------------------------------ 323 : * -- CALCUL DES VECTEURS DE LANCZOS -- 324 : * 324 : call gibtem(xkt) 325 : xkt1 = xkt1 + xkt 326 : * on cherche a obtenir Nmod modes Z sans depasser nvecma vecteurs X 326 : CALL LANCZO(IPKW2M,IPMASS,Nbmode,Nmod,Nvecma,Nmodopt 327 : c $ ,xktal1,iflu,nbonve,IPVALP,IPHI,IPERM,icyc1) 328 : $ ,xktal1,iflu,nbonve,IPVALP,IPHI,IPERM,icyc1,Northo) 329 : * on a obtenu nbonve candidats a devenir modes : IPVALP, X*PHI 330 : xkt1 = xkt1 + xktal1 331 : xktal1 = xktale 332 : IF(IERR.NE.0) RETURN 333 : 334 : * 335 : * ---------------------------------------------------- 336 : * -- W2FREQ POUR LES VALEURS PROPRES -- 337 : * -- RECOMBINAISON DES VECTEURS PROPRES Z = [X].PHI -- 338 : * -- CALCUL DU RESIDU (et test de convergence associée) -- 339 : * -- eventuellement, raffinement, purge des Z 340 : * et elimination des modes deja trouves -- 341 : * -- CREATION DE IPSOLU de ce cycle -- 342 : * 343 : Nramax = Northo 343 : CALL SIMUL3 (IPKW2M,IPMASS,IPVALP,IPHI,IPERM,Nmod,nbonve, 344 : $ Nramax,IPSOLU,LIMAG2) 345 : * on a retenu pour ce cycle nbonve modes 346 : * soit pour ce shift nbonZ modes stockés dans IPSOLU 347 : * et au total l ensemble des vecteurs modaux est stocké dans IPVECZ 348 : nbonZ = nbonZ + nbonve 348 : call gibtem(xkt3) 349 : xkt1 = xkt1 + xkt3 350 : if (iimpi.ge.6) then 351 : write(6,*) ' temps passé dans simul3 ' , xkt3 352 : write(6,*) ' Modes OK pour ce cycle, ce shift',nbonve,nbonZ 353 : $ ,' / ',Nmod 354 : endif 355 : IF(IERR.NE.0) RETURN 356 : 357 : 358 : C 359 : C --------------------------------------------------------- 360 : C -- CALCULS MASSE GEN. ET DEPL. GEN. ET NUMERO DE MODES -- 361 : c 362 : c (seulement si on a trouvé des modes supplémentaires...) 363 : if (nbonve.gt.0) then 363 : CALL SIMUL7(IPSOLU,IPMASS,IPKW2M,W2SHIF,IPLMOX,IPLMOY,IFLU) 363 : call gibtem(xkt7) 364 : xkt1 = xkt1 + xkt7 365 : if(iimpi.ge.6) write(6,*) ' temps passé dans simul7 ' , xkt7 366 : IF(IERR.NE.0) RETURN 367 : endif 368 : 369 : c petit message 370 : if (iimpi.ge.4) then 371 : WRITE(6,605) icycle,Nmod,frshif,nbonZ,(nvecZ + nbonZ) 372 : 605 FORMAT(2X,I5,2X,'|',2X,I5,8X,F12.3,3X,I5,10X,I5) 373 : endif 374 : 375 : c 376 : C --------------------------------------------------------- 377 : c -- faut-t-il cycler ? quel shift appliquer? -- 378 : c 379 : c --on a trouvé tous les modes necessaires pour ce shift : 380 : c on recommence un autre cycle avec un shift determine par strate 381 : c et pour Nmod modes (estime optimum % temps cpu) < Nmodma 382 : Nmod = min(Nmodopt,Nvecma/2) 383 : Nmodma = int(1.05*Nmod)+1 384 : c 385 : c -- strate : 386 : c -ajoute les modes de ce cycle (ipsolu) aux precedents (idist), 387 : c +si besoin, 388 : c -definit la direction de recherche de modes manquants, 389 : c -calcule un shift a appliquer au prochain cycle, 390 : c -verifie la completude du spectre 391 : c 392 : frtmp = frshif 393 : c 393 : call strate(IPSOLU,Nmodma,Nbmode,frtmp,Nmod,ifini 394 : $ ,shifti,idist,nbonZ,icycle,Northo 395 : $ ,IPRIGI,IPMASS,IPKW2M,nvp0,IFLU,incvrai) 396 : IF(IERR.NE.0) RETURN 397 : c 398 : * on a trouvé nbonZ modes pour ce shift qui s'ajoutent aux autres 399 : nvecZ = nvecZ + nbonZ 400 : 401 : c msolut = ipsolu 402 : * si strate decide qu on shifte... 403 : c if (frtmp.ne.frshif) then 404 : if ((frtmp.ne.frshif).and.(ifini.ne.1)) then 405 : frshif = frtmp 406 : W2SHIF= deuxpi * frshif 407 : W2SHIF= W2SHIF * W2SHIF 408 : lshift = .true. 409 : c il est possible que strate ait modifié tout seul IPKW2M 410 : c dans ce cas pas besoin de re-shifter (evite une factorisation) 411 : if (IPKWM0.ne.IPKW2M) then 412 : IPKWM0= IPKW2M 413 : cbp10/2010 : ipsolu mis a 0 a tous les cycles (pas seulement si on shifte) 414 : c mais nbonZ = nbre de mode pour ce shift (=plusieurs cycles) 415 : nbonZ = 0 416 : c IPSOLU = 0 417 : lshift = .false. 418 : endif 419 : endif 420 : 420 : call gibtem(xktstr) 421 : xkt1 = xkt1 + xktstr 422 : if(iimpi.ge.6) write(6,*) ' temps passé dans strate ' , xktstr 423 : 424 : c 425 : C --------------------------------------------------- 426 : c -- UN PEU DE MENAGE et SORTIE ? -- 427 : * 428 : * -- DESTRUCTION DES OBJETS DE TRAVAIL de ce shift -- 429 : * 430 : * on detruit les rigidites de ce cycle 430 : if(lshift.or.(ifini.eq.1)) CALL DTRIGI(IPKW2M) 431 : MLMOTS=IPLMOX 432 : MLMOT1=IPLMOY 433 : SEGSUP,MLMOTS,MLMOT1 434 : * 435 : * -- DESTRUCTION DES OBJETS DE TRAVAIL de ce cycle -- 436 : * 437 : * on detruit les vecteurs de Lanczos X de ce cycle 438 : mlent1=IPVECX 439 : segact mlent1 440 : do iou=1,mlent1.lect(/1) 441 : mvect1=mlent1.lect(iou) 442 : segsup,mvect1 443 : enddo 444 : segsup mlent1 445 : * 446 : * ---ON PENSE AVOIR GAGNÉ : ON QUITTE --- 447 : * 448 : if(ifini.eq.1) goto 999 449 : 450 : 451 : enddo 452 : *-----fin de la boucle des Cycles de Lanczos--------------------------* 453 : 454 : 455 : 999 continue 456 : * 457 : * -- DESTRUCTION DES OBJETS DE TRAVAIL communs a tous les cycles -- 458 : * 459 : * on detruit les modes convergés Z contenus dans IPVECZ 460 : mlent1=IPVECZ 461 : segact,mlent1 462 : do iou=1,mlent1.lect(/1) 463 : mvect1=mlent1.lect(iou) 464 : segsup,mvect1 465 : enddo 466 : segsup,mlent1 467 : IPVECZ = 0 468 : * on detruit les produits convergés M*Z/ZTMT contenus dans IPMZ 469 : mlent1=IPMZ 470 : segact,mlent1 471 : do iou=1,mlent1.lect(/1) 472 : mvect1=mlent1.lect(iou) 473 : segsup,mvect1 474 : enddo 475 : segsup,mlent1 476 : 476 : call gibtem(xkt) 477 : xkt1 = xkt1 + xkt 478 : write(6,*) ' temps passe dans VIBR SIMUL =',xkt1,' cs' 479 : if(iimpi.ge.4) 480 : & write(6,*) '=============================================' 481 : 482 : return 483 : * 484 : END 485 : 486 : 487 : 488 : 489 : 490 : 491 : 492 : 493 : 494 : 495 :
© Cast3M 2003 - All rights reserved.
Disclaimer