1 : C PTORI SOURCE CHAT 06/06/01 21:19:50 5450 2 : C CE SOUS-PROGRAMME RAMENE UN TORE SUR SES COORDONNEES PROPRES 3 : C LES DONNEES DU TORE SONT : LE CENTRE UN POINT DE L'AXE ET UN CENTRE 4 : C SECONDAIRE 5 : C 6 : SUBROUTINE PTORI(IOP,FER,XPROJ,NDEB,NUMNP,IPC,IPA,IP1,tcval,isens) 7 : IMPLICIT INTEGER(I-N) 8 : IMPLICIT REAL*8 (A-H,O-Z) 9 : -INC CCREEL 10 : -INC SMCOORD 11 : -INC CCOPTIO 12 : real*8 tcval(*) 13 : SEGMENT/FER/(NFI(ITT),MAI(IPP),ITOUR) 14 : SEGMENT XPROJ(3,IMAX) 15 : * tcval(1) 2 3 4 5 6 16 : * SAVE XORIG,YORIG,ZORIG,XAXE,YAXE,ZAXE 17 : * 7 8 9 18 : * SAVE XPC,YPC,ZPC 19 : * 10 11 20 : * SAVE ISENS,GRAY,PRAY 21 : IF (IOP.EQ.2) GOTO 100 22 : IMCT=MAI(ITOUR+1) 23 : INCT=MAI(1)+1 24 : IMAX=(IMCT**2)+10 24 : CALL LIRENT(IMAX,0,IRETOU) 25 : IF (IRETOU.NE.0) IMAX=MAX(1,IMAX) 26 : NDEB=IMCT+1 27 : SEGINI XPROJ 28 : SEGACT MCOORD 29 : C CENTRE DU TORE 30 : IREF=IPC*4-3 31 : XPC=XCOOR(IREF) 32 : YPC=XCOOR(IREF+1) 33 : ZPC=XCOOR(IREF+2) 34 : tcval(7)=xpc 35 : tcval(8)=ypc 36 : tcval(9)=zpc 37 : C POINT DE L'AXE 38 : IREF=IPA*4-3 39 : XPA=XCOOR(IREF) 40 : YPA=XCOOR(IREF+1) 41 : ZPA=XCOOR(IREF+2) 42 : C CENTRE SECONDAIRE 43 : IREF=IP1*4-3 44 : XP1=XCOOR(IREF) 45 : YP1=XCOOR(IREF+1) 46 : ZP1=XCOOR(IREF+2) 47 : XAXE=XPC-XPA 48 : YAXE=YPC-YPA 49 : ZAXE=ZPC-ZPA 50 : DAXE=SQRT(XAXE**2+YAXE**2+ZAXE**2) 50 : IF (DAXE.EQ.0.D0) CALL ERREUR(21) 51 : XAXE=XAXE/DAXE 52 : YAXE=YAXE/DAXE 53 : ZAXE=ZAXE/DAXE 54 : tcval(4)=xaxe 55 : tcval(5)=yaxe 56 : tcval(6)=zaxe 57 : C CALCUL DU GRAND RAYON 58 : GRAY=SQRT((XPC-XP1)**2+(YPC-YP1)**2+(ZPC-ZP1)**2) 59 : tcval(10)=gray 60 : IF (IIMPI.NE.0) WRITE (IOIMP,9901) GRAY 61 : 9901 FORMAT (/,' GRAND RAYON ',G12.5) 61 : IF (GRAY.EQ.0.D0) CALL ERREUR(21) 62 : C DEROULONS LE TORE (DANS LES DEUX SENS) 63 : IREF=4*NFI(IMCT)-3 64 : XV1=XCOOR(IREF)-XPC 65 : YV1=XCOOR(IREF+1)-YPC 66 : ZV1=XCOOR(IREF+2)-ZPC 67 : XI=XV1*XAXE+YV1*YAXE+ZV1*ZAXE 68 : XC1=XV1-XI*XAXE 69 : YC1=YV1-XI*YAXE 70 : ZC1=ZV1-XI*ZAXE 71 : C CENTRE SECONDAIRE 72 : SC1=SQRT(XC1**2+YC1**2+ZC1**2) 73 : U1=XC1/SC1 74 : V1=YC1/SC1 75 : W1=ZC1/SC1 76 : XD1=XPC+U1*GRAY 77 : YD1=YPC+V1*GRAY 78 : ZD1=ZPC+W1*GRAY 79 : XORIG=U1 80 : YORIG=V1 81 : ZORIG=W1 82 : tcval(1)=xorig 83 : tcval(2)=yorig 84 : tcval(3)=zorig 85 : U2=YAXE*W1-ZAXE*V1 86 : V2=ZAXE*U1-XAXE*W1 87 : W2=XAXE*V1-YAXE*U1 88 : XVOR1=XV1+XPC-XD1 89 : YVOR1=YV1+YPC-YD1 90 : ZVOR1=ZV1+ZPC-ZD1 91 : COT2=ATAN2(U2*(YAXE*ZVOR1-ZAXE*YVOR1)+V2* 92 : # (ZAXE*XVOR1-XAXE*ZVOR1)+W2*(XAXE*YVOR1-YAXE*XVOR1), 93 : # XAXE*XVOR1+YAXE*YVOR1+ZAXE*ZVOR1) 94 : RAYON=0.D0 95 : COT1=0.D0 96 : DO 1 I=INCT,max(IMCT,mai(itour+2)) 97 : II=NFI(I) 98 : IREF=4*II-3 99 : XV2=XCOOR(IREF)-XPC 100 : YV2=XCOOR(IREF+1)-YPC 101 : ZV2=XCOOR(IREF+2)-ZPC 102 : XI=XV2*XAXE+YV2*YAXE+ZV2*ZAXE 103 : XC2=XV2-XI*XAXE 104 : YC2=YV2-XI*YAXE 105 : ZC2=ZV2-XI*ZAXE 106 : ANG=ATAN2(XAXE*(YC1*ZC2-YC2*ZC1)+YAXE*(ZC1*XC2-ZC2*XC1)+ 107 : # ZAXE*(XC1*YC2-XC2*YC1),XC1*XC2+YC1*YC2+ZC1*ZC2) 108 : COT1=COT1+ANG 109 : XPROJ(1,I)=COT1*GRAY 110 : C CENTRE SECONDAIRE 111 : SC2=SQRT(XC2**2+YC2**2+ZC2**2) 111 : if (sc2.eq.0.D0) call erreur(21) 112 : U1=XC2/SC2 113 : V1=YC2/SC2 114 : W1=ZC2/SC2 115 : XD2=XPC+U1*GRAY 116 : YD2=YPC+V1*GRAY 117 : ZD2=ZPC+W1*GRAY 118 : U2=YAXE*W1-ZAXE*V1 119 : V2=ZAXE*U1-XAXE*W1 120 : W2=XAXE*V1-YAXE*U1 121 : XVOR2=XV2+XPC-XD2 122 : YVOR2=YV2+YPC-YD2 123 : ZVOR2=ZV2+ZPC-ZD2 124 : ANG=ATAN2(U2*(YAXE*ZVOR2-ZAXE*YVOR2)+V2* 125 : # (ZAXE*XVOR2-XAXE*ZVOR2)+W2*(XAXE*YVOR2-YAXE*XVOR2), 126 : # XAXE*XVOR2+YAXE*YVOR2+ZAXE*ZVOR2) 127 : ADIF=ANG-COT2 128 : ADIF=ADIF-INT(ADIF/(2*XPI))*2*XPI 129 : IF (ADIF.GE.XPI) ADIF=ADIF-2*XPI 130 : IF (ADIF.LE.-XPI) ADIF=ADIF+2*XPI 131 : COT2=COT2+ADIF 132 : XPROJ(2,I)=COT2 133 : XPROJ(3,I)=XCOOR(IREF+3) 134 : if (i.le.imct) RAYON=RAYON+XVOR2**2+YVOR2**2+ZVOR2**2 135 : XC1=XC2 136 : YC1=YC2 137 : ZC1=ZC2 138 : 1 CONTINUE 139 : RAYON2=RAYON/(IMCT-INCT+1) 140 : PRAY=SQRT(RAYON2) 141 : tcval(11)=pray 142 : IF (IIMPI.NE.0) WRITE (IOIMP,9902) PRAY 143 : 9902 FORMAT (' PETIT RAYON ',G12.5) 144 : DO 2 I=INCT,max(IMCT,mai(itour+2)) 145 : XPROJ(2,I)=XPROJ(2,I)*PRAY 146 : II=NFI(I) 147 : NFI(I)=I 148 : IREF=4*II-3 149 : XV=XCOOR(IREF)-XPC 150 : YV=XCOOR(IREF+1)-YPC 151 : ZV=XCOOR(IREF+2)-ZPC 152 : XI=XV*XAXE+YV*YAXE+ZV*ZAXE 153 : XC=XV-XI*XAXE 154 : YC=YV-XI*YAXE 155 : ZC=ZV-XI*ZAXE 156 : C CENTRE SECONDAIRE 157 : SC=SQRT(XC**2+YC**2+ZC**2) 158 : U1=XC/SC 159 : V1=YC/SC 160 : W1=ZC/SC 161 : XC=XPC+U1*GRAY 162 : YC=YPC+V1*GRAY 163 : ZC=ZPC+W1*GRAY 164 : XVOR=XV+XPC-XC 165 : YVOR=YV+YPC-YC 166 : ZVOR=ZV+ZPC-ZC 167 : RAY2=XVOR**2+YVOR**2+ZVOR**2 168 : RAP=RAY2/RAYON2 169 : IF (i.gt.imct.or.(RAP.LE.1.02D0.AND.RAP.GE.0.98D0)) GOTO 2 170 : IF (IIMPI.NE.0) WRITE (IOIMP,9903) I,RAP 171 : 9903 FORMAT(' POINT ',I6,' ERREUR RELATIVE SUR LE PETIT RAYON ',G12 172 : $ .5) 172 : CALL ERREUR(21) 173 : RETURN 174 : 2 CONTINUE 175 : C IL FAUT TOURNER DANS LE BON SENS 176 : SURF=0.D0 177 : DO 3 IT=1,ITOUR 178 : II1=MAI(IT-1+1)+1 179 : II2=MAI(IT+1) 180 : XV1=XPROJ(1,II2) 181 : YV1=XPROJ(2,II2) 182 : DO 31 I=II1,II2 183 : XV2=XPROJ(1,I) 184 : YV2=XPROJ(2,I) 185 : IF (XPROJ(3,I).EQ.0.D0) XPROJ(3,I)=SQRT((XV2-XV1)**2+(YV2 186 : $ -YV1)**2) 187 : SURF=SURF+XV1*YV2-XV2*YV1 188 : XV1=XV2 189 : YV1=YV2 190 : 31 CONTINUE 191 : 3 CONTINUE 192 : ISENS=1 193 : IF (SURF.GT.0.D0) GOTO 5 194 : ISENS=-1 195 : DO 4 I=INCT,max(IMCT,mai(itour+2)) 196 : XPROJ(1,I)=-XPROJ(1,I) 197 : 4 CONTINUE 198 : 5 CONTINUE 199 : RETURN 200 : C TRANSFORMATION INVERSE 201 : 100 CONTINUE 202 : xorig=tcval(1) 203 : yorig=tcval(2) 204 : zorig=tcval(3) 205 : xaxe=tcval(4) 206 : yaxe=tcval(5) 207 : zaxe=tcval(6) 208 : xpc=tcval(7) 209 : ypc=tcval(8) 210 : zpc=tcval(9) 211 : gray= tcval(10) 212 : pray=tcval(11) 213 : SEGACT MCOORD 214 : XREP1=XORIG 215 : YREP1=YORIG 216 : ZREP1=ZORIG 217 : XREP2=YAXE*ZREP1-ZAXE*YREP1 218 : YREP2=ZAXE*XREP1-XAXE*ZREP1 219 : ZREP2=XAXE*YREP1-YAXE*XREP1 220 : IF (NDEB.GT.NUMNP) GOTO 102 221 : NBPTA=XCOOR(/1)/(IDIM+1) 222 : NBPTS=NBPTA+NUMNP-NDEB+1 223 : SEGADJ MCOORD 224 : DO 101 I=NDEB,NUMNP 225 : ANG=XPROJ(1,I)*ISENS/GRAY 226 : SI=SIN(ANG) 227 : CO=COS(ANG) 228 : XC1=(XREP1*CO+XREP2*SI)*GRAY 229 : YC1=(YREP1*CO+YREP2*SI)*GRAY 230 : ZC1=(ZREP1*CO+ZREP2*SI)*GRAY 231 : SC1=SQRT(XC1**2+YC1**2+ZC1**2) 232 : ANG=XPROJ(2,I)/PRAY 233 : SI=SIN(ANG) 234 : CO=COS(ANG) 235 : XV1=(XAXE*CO+XC1*SI/SC1)*PRAY 236 : YV1=(YAXE*CO+YC1*SI/SC1)*PRAY 237 : ZV1=(ZAXE*CO+ZC1*SI/SC1)*PRAY 238 : XCOOR(NBPTA*(IDIM+1)+1)=XPC+XC1+XV1 239 : XCOOR(NBPTA*(IDIM+1)+2)=YPC+YC1+YV1 240 : XCOOR(NBPTA*(IDIM+1)+3)=ZPC+ZC1+ZV1 241 : XCOOR((NBPTA+1)*(IDIM+1))=XPROJ(3,I) 242 : NBPTA=NBPTA+1 243 : 101 CONTINUE 244 : 102 CONTINUE 245 : SEGSUP XPROJ 246 : RETURN 247 : END 248 : 249 : 250 : 251 : 252 : 253 :
© Cast3M 2003 - All rights reserved.
Disclaimer