1 : C PSPHE SOURCE CHAT 06/06/01 21:19:43 5450 2 : C CE SOUS-PROGRAMME RAMENNE UNE SPHERE SUR SES COORDONNEES PROPRES 3 : C 4 : SUBROUTINE PSPHE(IOP,FER,XPROJ,NDEB,NUMNP,ICENT,tcval) 5 : IMPLICIT INTEGER(I-N) 6 : IMPLICIT REAL*8 (A-H,O-Z) 7 : -INC SMCOORD 8 : -INC CCOPTIO 9 : real*8 tcval(*) 10 : SEGMENT/FER/(NFI(ITT),MAI(IPP),ITOUR) 11 : SEGMENT XPROJ(3,IMAX) 12 : SEGMENT XINT(3,max(MAI(ITOUR+1),mai(itour+2))) 13 : * tcval (1) 2 3 4 5 6 7 8 9 14 : * SAVE XVEC1,YVEC1,ZVEC1,XVEC2,YVEC2,ZVEC2,XGRAV,YGRAV,ZGRAV 15 : * tcval (10) 11 12 13 16 : * SAVE RAYNV,XINV,YINV,ZINV 17 : IF (IOP.EQ.2) GOTO 100 18 : IMCT=MAI(ITOUR+1) 19 : INCT=MAI(1)+1 20 : IMAX=(IMCT**2)/4+10 20 : CALL LIRENT(IMAX,0,IRETOU) 21 : IF (IRETOU.NE.0) IMAX=MAX(1,IMAX) 22 : NDEB=IMCT+1 23 : SEGINI XPROJ 24 : SEGINI XINT 25 : SEGACT MCOORD 26 : C CENTRE DE LA SPHERE 27 : IREF=ICENT*4-3 28 : XCENT=XCOOR(IREF) 29 : YCENT=XCOOR(IREF+1) 30 : ZCENT=XCOOR(IREF+2) 31 : C CENTRE DE GRAVITE POUR DETERMINER LE CENTRE DE L'INVERSION 32 : C CALCUL DU RAYON DE LA SPHERE 33 : XGRAVS=0 34 : YGRAVS=0 35 : ZGRAVS=0 36 : RAYON=0. 37 : DO 1 I=INCT,IMCT 38 : IREF=NFI(I)*4-3 39 : XP=XCOOR(IREF) 40 : YP=XCOOR(IREF+1) 41 : ZP=XCOOR(IREF+2) 42 : XGRAVS=XGRAVS+XP 43 : YGRAVS=YGRAVS+YP 44 : ZGRAVS=ZGRAVS+ZP 45 : RAYON=RAYON+(XCENT-XP)**2+(YCENT-YP)**2+(ZCENT-ZP)**2 46 : 1 CONTINUE 47 : XGRAVS=XGRAVS/(IMCT-INCT+1) 48 : YGRAVS=YGRAVS/(IMCT-INCT+1) 49 : ZGRAVS=ZGRAVS/(IMCT-INCT+1) 50 : RAYON=SQRT(RAYON/(IMCT-INCT+1)) 51 : 52 : C CENTRE DE L'INVERSION 53 : XDIR=XCENT-XGRAVS 54 : YDIR=YCENT-YGRAVS 55 : ZDIR=ZCENT-ZGRAVS 56 : DDIR=SQRT(XDIR**2+YDIR**2+ZDIR**2) 57 : COF=RAYON/DDIR 58 : XINV=XCENT+COF*XDIR 59 : YINV=YCENT+COF*YDIR 60 : ZINV=ZCENT+COF*ZDIR 61 : tcval(11)=xinv 62 : tcval(12)=yinv 63 : tcval(13)=zinv 64 : C EN AVANT POUR L'INVERSION ATTENTION AU CALCUL DE LA DENSITE 65 : RAYNV=4*RAYON**2 66 : tcval(10)=raynv 67 : DO 40 I=INCT,max(IMCT,mai(itour+2)) 68 : II=NFI(I) 69 : IREF=II*4-3 70 : XV=XCOOR(IREF)-XINV 71 : YV=XCOOR(IREF+1)-YINV 72 : ZV=XCOOR(IREF+2)-ZINV 73 : DV=XV**2+YV**2+ZV**2 73 : IF (DV.EQ.0.) CALL ERREUR(21) 74 : IF (IERR.NE.0) RETURN 75 : XINT(1,I)=XV*RAYNV/DV 76 : XINT(2,I)=YV*RAYNV/DV 77 : XINT(3,I)=ZV*RAYNV/DV 78 : 40 CONTINUE 79 : C CENTRE DE GRAVITE DU PLAN DE PROJECTION 80 : XGRAV=0 81 : YGRAV=0 82 : ZGRAV=0 83 : DO 51 I=INCT,IMCT 84 : XGRAV=XGRAV+XINT(1,I) 85 : YGRAV=YGRAV+XINT(2,I) 86 : ZGRAV=ZGRAV+XINT(3,I) 87 : 51 CONTINUE 88 : XGRAV=XGRAV/(IMCT-INCT+1) 89 : YGRAV=YGRAV/(IMCT-INCT+1) 90 : ZGRAV=ZGRAV/(IMCT-INCT+1) 91 : tcval(7)=xgrav 92 : tcval(8)=ygrav 93 : tcval(9)=zgrav 94 : C VECTEUR NORMAL 95 : XNORM=0 96 : YNORM=0 97 : ZNORM=0 98 : DO 52 IT=1,ITOUR 99 : IPR=MAI(IT+1) 100 : XV1=XINT(1,IPR)-XGRAV 101 : YV1=XINT(2,IPR)-YGRAV 102 : ZV1=XINT(3,IPR)-ZGRAV 103 : DO 52 I=MAI(IT-1+1)+1,MAI(IT+1) 104 : XV2=XINT(1,I)-XGRAV 105 : YV2=XINT(2,I)-YGRAV 106 : ZV2=XINT(3,I)-ZGRAV 107 : XNORM=XNORM+YV1*ZV2-ZV1*YV2 108 : YNORM=YNORM+ZV1*XV2-XV1*ZV2 109 : ZNORM=ZNORM+XV1*YV2-XV2*YV1 110 : XV1=XV2 111 : YV1=YV2 112 : ZV1=ZV2 113 : 52 CONTINUE 114 : DNORM=SQRT (XNORM**2+YNORM**2+ZNORM**2) 115 : XNORM=XNORM/DNORM 116 : YNORM=YNORM/DNORM 117 : ZNORM=ZNORM/DNORM 118 : C FORMATION DU REPERE 119 : IF (ABS(XNORM).LT.0.1) THEN 120 : XVEC1=1-XNORM*XNORM 121 : YVEC1= -XNORM*YNORM 122 : ZVEC1= -XNORM*ZNORM 123 : ELSE 124 : XVEC1= -YNORM*XNORM 125 : YVEC1=1-YNORM*YNORM 126 : ZVEC1= -YNORM*ZNORM 127 : ENDIF 128 : DVEC1=XVEC1**2+YVEC1**2+ZVEC1**2 129 : DVEC1=SQRT(DVEC1) 130 : XVEC1=XVEC1/DVEC1 131 : YVEC1=YVEC1/DVEC1 132 : ZVEC1=ZVEC1/DVEC1 133 : XVEC2=YNORM*ZVEC1-ZNORM*YVEC1 134 : YVEC2=ZNORM*XVEC1-XNORM*ZVEC1 135 : ZVEC2=XNORM*YVEC1-YNORM*XVEC1 136 : tcval(1)=xvec1 137 : tcval(2)=yvec1 138 : tcval(3)=Zvec1 139 : tcval(4)=xvec2 140 : tcval(5)=Yvec2 141 : tcval(6)=Zvec2 142 : C EN AVANT POUR LA PROJECTION 143 : DO 53 I=INCT,max(IMCT,mai(itour+2)) 144 : II=NFI(I) 145 : NFI(I)=I 146 : XRE=XINT(1,I)-XGRAV 147 : YRE=XINT(2,I)-YGRAV 148 : ZRE=XINT(3,I)-ZGRAV 149 : XPROJ(1,I)=XRE*XVEC1+YRE*YVEC1+ZRE*ZVEC1 150 : XPROJ(2,I)=XRE*XVEC2+YRE*YVEC2+ZRE*ZVEC2 151 : XTEST =XRE*XNORM+YRE*YNORM+ZRE*ZNORM 151 : IF (ABS(XTEST).GT.DNORM*1E-2.and.i.le.imct) CALL ERREUR(21) 152 : IF (IERR.NE.0) RETURN 153 : XPROJ(3,I)=XCOOR(II*4) 154 : 53 CONTINUE 155 : C VERIFIER LES DENSITES 156 : DO 54 IT=1,ITOUR 157 : II1=MAI(IT-1+1)+1 158 : II2=MAI(IT+1) 159 : DO 54 I=II1,II2 160 : IF (XPROJ(3,I).NE.0) GOTO 54 161 : IAP=I+1 162 : IF (IAP.GT.II2) IAP=II1 163 : XPROJ(3,I)=SQRT((XPROJ(1,I)-XPROJ(1,IAP))**2+(XPROJ(2,I)-XPROJ(2, 164 : # IAP))**2) 165 : 54 CONTINUE 166 : SEGSUP XINT 167 : RETURN 168 : 100 CONTINUE 169 : C ON RECONSTITUE LE MAILLAGE 170 : xvec1=tcval(1) 171 : yvec1=tcval(2) 172 : zvec1=tcval(3) 173 : xvec2=tcval(4) 174 : yvec2=tcval(5) 175 : zvec2=tcval(6) 176 : xgrav=tcval(7) 177 : ygrav=tcval(8) 178 : zgrav=tcval(9) 179 : raynv=tcval(10) 180 : xinv=tcval(11) 181 : yinv=tcval(12) 182 : zinv=tcval(13) 183 : SEGACT MCOORD 184 : IF (NDEB.GT.NUMNP) GOTO 111 185 : NBPTA=XCOOR(/1)/(IDIM+1) 186 : NBPTS=NBPTA+NUMNP-NDEB+1 187 : SEGADJ MCOORD 188 : DO 110 I=NDEB,NUMNP 189 : XI=XPROJ(1,I)*XVEC1+XPROJ(2,I)*XVEC2+XGRAV 190 : YI=XPROJ(1,I)*YVEC1+XPROJ(2,I)*YVEC2+YGRAV 191 : ZI=XPROJ(1,I)*ZVEC1+XPROJ(2,I)*ZVEC2+ZGRAV 192 : DI=XI**2+YI**2+ZI**2 193 : XCOOR(NBPTA*(IDIM+1)+1)=XI*RAYNV/DI+XINV 194 : XCOOR(NBPTA*(IDIM+1)+2)=YI*RAYNV/DI+YINV 195 : XCOOR(NBPTA*(IDIM+1)+3)=ZI*RAYNV/DI+ZINV 196 : XCOOR((NBPTA+1)*(IDIM+1))=XPROJ(3,I) 197 : NBPTA=NBPTA+1 198 : 110 CONTINUE 199 : 111 CONTINUE 200 : SEGSUP XPROJ 201 : RETURN 202 : END 203 : 204 : 205 :
© Cast3M 2003 - All rights reserved.
Disclaimer