Télécharger ptori.eso

Retour à la liste

Numérotation des lignes :

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

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