Télécharger ptori.eso

Retour à la liste

Numérotation des lignes :

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales