1 : C PCONE SOURCE CHAT 06/06/01 21:19:10 5450 2 : C CE SOUS-PROGRAMME RAMENNE UN CONE SUR SES COORDONNEES PROPRES 3 : C 4 : SUBROUTINE PCONE(IOP,FER,XPROJ,NDEB,NUMNP,ISOM,IP2,tcval,isens) 5 : IMPLICIT INTEGER(I-N) 6 : IMPLICIT REAL*8 (A-H,O-Z) 7 : -INC SMCOORD 8 : -INC CCOPTIO 9 : SEGMENT/FER/(NFI(ITT),MAI(IPP),ITOUR) 10 : SEGMENT XPROJ(3,IMAX) 11 : real*8 tcval(*) 12 : * tcval(1) 2 3 4 5 6 7 13 : * SAVE XORIG,YORIG,ZORIG,XAXE,YAXE,ZAXE,ISENS,SANGS 14 : * tcval(8) 9 10 11 15 : * SAVE XSOM,YSOM,ZSOM,CANGS 16 : IF (IOP.EQ.2) GOTO 100 17 : IMCT=MAI(ITOUR+1) 18 : INCT=MAI(1)+1 19 : IMAX=IMCT**2 19 : CALL LIRENT(IMAX,0,IRETOU) 20 : IF (IRETOU.NE.0) IMAX=MAX(1,IMAX) 21 : NDEB=IMCT+1 22 : SEGINI XPROJ 23 : SEGACT MCOORD 24 : C AXE DU CONE 25 : IREF=ISOM*4-3 26 : XSOM=XCOOR(IREF) 27 : YSOM=XCOOR(IREF+1) 28 : ZSOM=XCOOR(IREF+2) 29 : tcval(8)=xsom 30 : tcval(9)=ysom 31 : tcval(10)=Zsom 32 : IREF=IP2*4-3 33 : XP2=XCOOR(IREF) 34 : YP2=XCOOR(IREF+1) 35 : ZP2=XCOOR(IREF+2) 36 : XAXE=XP2-XSOM 37 : YAXE=YP2-YSOM 38 : ZAXE=ZP2-ZSOM 39 : DAXE=SQRT(XAXE**2+YAXE**2+ZAXE**2) 39 : IF (DAXE.EQ.0.) CALL ERREUR(21) 40 : IF (IERR.NE.0) RETURN 41 : XAXE=XAXE/DAXE 42 : YAXE=YAXE/DAXE 43 : ZAXE=ZAXE/DAXE 44 : tcval(4)=xaxe 45 : tcval(5)=yaxe 46 : tcval(6)=zaxe 47 : C DEROULONS LE CONE 48 : IREF=4*NFI(IMCT)-3 49 : XV1=XCOOR(IREF)-XSOM 50 : YV1=XCOOR(IREF+1)-YSOM 51 : ZV1=XCOOR(IREF+2)-ZSOM 52 : PV1=XV1*XAXE+YV1*YAXE+ZV1*ZAXE 53 : XV1=XV1-PV1*XAXE 54 : YV1=YV1-PV1*YAXE 55 : ZV1=ZV1-PV1*ZAXE 56 : DV1=SQRT(XV1**2+YV1**2+ZV1**2) 56 : IF (DV1.EQ.0.) CALL ERREUR(21) 57 : IF (IERR.NE.0) RETURN 58 : XV1=XV1/DV1 59 : YV1=YV1/DV1 60 : ZV1=ZV1/DV1 61 : XORIG=XV1 62 : YORIG=YV1 63 : ZORIG=ZV1 64 : tcval(1)=xorig 65 : tcval(2)=yorig 66 : tcval(3)=zorig 67 : ANGSOM=0 68 : COT=0 69 : DO 1 I=INCT,max(IMCT,mai(itour+2)) 70 : II=NFI(I) 71 : IREF=4*II-3 72 : XV2=XCOOR(IREF)-XSOM 73 : YV2=XCOOR(IREF+1)-YSOM 74 : ZV2=XCOOR(IREF+2)-ZSOM 75 : XPROJ(1,I)=SQRT(XV2**2+YV2**2+ZV2**2) 76 : XPROJ(3,I)=XCOOR(IREF+3) 77 : XI=XPROJ(1,I) 78 : C XI EST LE RAYON EN COORDONNEES POLAIRES QU'EST L'ANGLE? 79 : PROJ=XV2*XAXE+YV2*YAXE+ZV2*ZAXE 80 : if (i.le.imct) ANGSOM=ANGSOM+PROJ/XI 81 : XV2=XV2-PROJ*XAXE 82 : YV2=YV2-PROJ*YAXE 83 : ZV2=ZV2-PROJ*ZAXE 84 : DV2=SQRT(XV2**2+YV2**2+ZV2**2) 84 : IF (i.le.imct.and.DV2.EQ.0.) CALL ERREUR(21) 85 : IF (IERR.NE.0) RETURN 86 : XV2=XV2/DV2 87 : YV2=YV2/DV2 88 : ZV2=ZV2/DV2 89 : ANG=ATAN2(XAXE*(YV1*ZV2-ZV1*YV2)+YAXE*(ZV1*XV2-XV1*ZV2)+ZAXE* 90 : # (XV1*YV2-YV1*XV2),XV1*XV2+YV1*YV2+ZV1*ZV2) 91 : COT=COT+ANG 92 : XPROJ(2,I)=COT 93 : XV1=XV2 94 : YV1=YV2 95 : ZV1=ZV2 96 : 1 CONTINUE 97 : CANGS=ANGSOM/(IMCT-INCT+1) 98 : SANGS=SQRT(1-CANGS**2) 99 : tcval(7)=sangs 100 : tcval(11)=cangs 101 : DO 2 I=INCT,max(IMCT,mai(itour+2)) 102 : XPROJ(2,I)=XPROJ(2,I)*SANGS 103 : RR=XPROJ(1,I) 104 : XPROJ(1,I)=RR*COS(XPROJ(2,I)) 105 : XPROJ(2,I)=RR*SIN(XPROJ(2,I)) 106 : II=NFI(I) 107 : NFI(I)=I 108 : IREF=4*II-3 109 : XV=XCOOR(IREF)-XSOM 110 : YV=XCOOR(IREF+1)-YSOM 111 : ZV=XCOOR(IREF+2)-ZSOM 112 : SCAL=XV*XAXE+YV*YAXE+ZV*ZAXE 113 : DV=SQRT(XV**2+YV**2+ZV**2) 114 : CAN=SCAL/DV 115 : RAP=CANGS/CAN 115 : IF (i.le.imct.and.(RAP.GT.1.01.OR.RAP.LT.0.99)) CALL ERREUR(21) 116 : 2 CONTINUE 117 : C IL FAUT TOURNER DANS LE BON SENS 118 : SURF=0 119 : DO 3 IT=1,ITOUR 120 : II1=MAI(IT-1+1)+1 121 : II2=MAI(IT+1) 122 : XV1=XPROJ(1,II2) 123 : YV1=XPROJ(2,II2) 124 : DO 3 I=II1,II2 125 : XV2=XPROJ(1,I) 126 : YV2=XPROJ(2,I) 127 : IF (XPROJ(3,I).EQ.0) XPROJ(3,I)=SQRT((XV2-XV1)**2+(YV2-YV1)**2) 128 : SURF=SURF+XV1*YV2-XV2*YV1 129 : XV1=XV2 130 : YV1=YV2 131 : 3 CONTINUE 132 : ISENS=1 133 : IF (SURF.GT.0.) GOTO 5 134 : ISENS=-1 135 : DO 4 I=INCT,max(IMCT,mai(itour+2)) 136 : XPROJ(1,I)=-XPROJ(1,I) 137 : 4 CONTINUE 138 : 5 CONTINUE 139 : RETURN 140 : C TRANSFORMATION INVERSE 141 : 100 CONTINUE 142 : xorig=tcval(1) 143 : yorig=tcval(2) 144 : zorig=tcval(3) 145 : xaxe=tcval(4) 146 : yaxe=tcval(5) 147 : zaxe=tcval(6) 148 : sangs=tcval(7) 149 : xsom=tcval(8) 150 : ysom=tcval(9) 151 : zsom=tcval(10) 152 : cangs=tcval(11) 153 : SEGACT MCOORD 154 : XREP1=XORIG 155 : YREP1=YORIG 156 : ZREP1=ZORIG 157 : XREP2=YAXE*ZREP1-ZAXE*YREP1 158 : YREP2=ZAXE*XREP1-XAXE*ZREP1 159 : ZREP2=XAXE*YREP1-YAXE*XREP1 160 : IF (NUMNP.LT.NDEB) GOTO 102 161 : NBPTA=XCOOR(/1)/4 162 : NBPTS=NBPTA+NUMNP-NDEB+1 163 : SEGADJ MCOORD 164 : DO 101 I=NDEB,NUMNP 165 : R=SQRT(XPROJ(1,I)**2+XPROJ(2,I)**2) 165 : IF (R.EQ.0.) CALL ERREUR(21) 166 : IF (IERR.NE.0.) RETURN 167 : ANG=ATAN2(XPROJ(2,I),XPROJ(1,I)*ISENS) 168 : XPROJ(1,I)=R 169 : XPROJ(2,I)=ANG/SANGS 170 : SI=SIN(XPROJ(2,I)) 171 : CO=COS(XPROJ(2,I)) 172 : XCOOR(NBPTA*4+1)=(XREP1*CO+XREP2*SI)*R*SANGS+R*CANGS*XAXE+XSOM 173 : XCOOR(NBPTA*4+2)=(YREP1*CO+YREP2*SI)*R*SANGS+R*CANGS*YAXE+YSOM 174 : XCOOR(NBPTA*4+3)=(ZREP1*CO+ZREP2*SI)*R*SANGS+R*CANGS*ZAXE+ZSOM 175 : XCOOR(NBPTA*4+4)=XPROJ(3,I) 176 : NBPTA=NBPTA+1 177 : 101 CONTINUE 178 : 102 CONTINUE 179 : SEGSUP XPROJ 180 : RETURN 181 : END 182 : 183 :
© Cast3M 2003 - All rights reserved.
Disclaimer