Télécharger inters.eso

Retour à la liste

Numérotation des lignes :

  1. C INTERS SOURCE BP208322 16/11/18 21:17:49 9177
  2. SUBROUTINE INTERS
  3. C-----------------------------------------------------------------------
  4. C
  5. C INTERSECTIONS DE SURFACES (PLAN,SPHERE,CYLINDRE,CONE,TORE)
  6. C
  7. C-----------------------------------------------------------------------
  8. C
  9. C IOB1 EVENTUELLEMENT OBJET QUE L'ON PROLONGE
  10. C IOB2 EVENTUELLEMENT OBJET QUE L'ON AVANCE
  11. C IPOIN1 POINT INITIAL
  12. C IPOIN2 POINT FINAL
  13. C ICLE TYPE DE LA PREMIERE SURFACE
  14. C JCLE TYPE DE LA DEUXIEME SURFACE
  15. C XPAUX TABLEAU DE POINTS AUXILIAIRES SERVANT A RECTIFIER L'INTERSE
  16. C LA QUATRIEME COMPOSANTE EST L'ABCISSE CURVILIGNE
  17. C
  18. C-----------------------------------------------------------------------
  19. C
  20. C 17/09/93 : intersection de deux MELEMEs (sens ensembliste)
  21. C
  22. C PM 30/08/05 : Erreur si intersection vide, sauf si mot-clef 'NOVERIF'
  23. C
  24. C BP 28/08/2012 : Ajout de l'intersection geometrique de 2 maillages
  25. C
  26. C-----------------------------------------------------------------------
  27. c
  28. IMPLICIT INTEGER(I-N)
  29. IMPLICIT REAL*8 (A-H,O-Z)
  30.  
  31. -INC CCOPTIO
  32. -INC CCGEOME
  33.  
  34. -INC SMCOORD
  35. -INC SMELEME
  36.  
  37. SEGMENT TABPAR(NPOIN)
  38.  
  39. LOGICAL VERIF,ltelq
  40. CHARACTER*4 MCLE(8),MCLE2(1),MCLE3(1)
  41. C INITIALISATION A L'USAGE DE L'OPTIMISEUR
  42. PARAMETER ( NBT = 22 )
  43. DIMENSION XPAUX(4,NBT)
  44.  
  45. DATA RAYI/1./
  46. DATA MCLE/'PLAN','SPHE','CYLI','CONI','TORI','DINI','DFIN','LONG'/
  47. DATA MCLE2/'NOVE'/
  48. DATA MCLE3/'GEOM'/
  49. C
  50. C-----------------------------------------------------------------------
  51. C INTERSECTION DE DEUX MAILLAGES
  52. C
  53. C sens ensembliste ou geometrique?
  54. CALL LIRMOT(MCLE3,1,IGEO3,0)
  55. c write(ioimp,*) 'option GEOM ? ',IGEO3
  56. C
  57. c lecture des maillages
  58. CALL LIROBJ('MAILLAGE',IOB1,0,IRETOU)
  59. IF (IRETOU.EQ.1) THEN
  60. CALL LIROBJ('MAILLAGE',IOB2,0,IRETOV)
  61. C On laisse la syntaxe NOVERIF pour l'intersection Géométrique
  62. CALL LIRMOT(MCLE2,1,ICLE,0)
  63. VERIF=(ICLE.EQ.0)
  64.  
  65. C------ Intersection ensembliste ------------------
  66. if(IGEO3.eq.0) then
  67. IF (IRETOV.EQ.1) THEN
  68. C Intersection de deux maillages
  69. CALL INTERB(IOB1,IOB2,IRETIB,IPT3)
  70. IF (IRETIB.EQ.1) THEN
  71. C On fournit un maillage vide
  72. NBNN=0
  73. NBELEM=0
  74. NBSOUS=0
  75. NBREF=0
  76. SEGINI,IPT3
  77. IPT3.ITYPEL=0
  78. SEGDES,IPT3
  79. ENDIF
  80. CALL ECROBJ('MAILLAGE',IPT3)
  81. RETURN
  82. ELSE
  83. CALL REFUS
  84. ENDIF
  85.  
  86. C------ Intersection geometrique ------------------
  87. else
  88. IOB3=0
  89. IOB4=0
  90. CALL INTGEO(IOB1,IOB2,IOB3,IOB4,ICLE)
  91. C Cas ou n'a pas pu calculer l'intersection
  92. if(IOB3.eq.0) then
  93. C On fournit 2 maillage vide + les 2 initiaux
  94. NBNN=0
  95. NBELEM=0
  96. NBSOUS=0
  97. NBREF=0
  98. SEGINI,IPT3
  99. IPT3.ITYPEL=0
  100. IOB3=IPT3
  101. SEGDES,IPT3
  102. IF(IRETOV.EQ.1)THEN
  103. SEGINI,IPT3,IPT4
  104. IPT4.ITYPEL=0
  105. IOB4=IPT4
  106. SEGDES,IPT4
  107. ENDIF
  108. endif
  109. IF(IRETOV.EQ.1) CALL ECROBJ('MAILLAGE',IOB4)
  110. CALL ECROBJ('MAILLAGE',IOB3)
  111. IF(IRETOV.EQ.1) CALL ECROBJ('MAILLAGE',IOB2)
  112. CALL ECROBJ('MAILLAGE',IOB1)
  113. RETURN
  114. endif
  115.  
  116. ENDIF
  117.  
  118.  
  119. C-----------------------------------------------------------------------
  120. C CREATION D'UNE LIGNE RESULTAT D'UNE INTERSECTION DE SURFACES
  121. C DEFINIES ANALYTIQUEMENT (PLAN, SPHERE, ...)
  122. C
  123. IF (IDIM.NE.3) CALL ERREUR(14)
  124. IF (IERR.NE.0) RETURN
  125.  
  126. C Y A T'IL UN FACTEUR DE DECOUPAGE
  127. INBR=0
  128. CALL LIRENT(INBR,0,IRETOU)
  129. IMPOI=0
  130. IMPOF=0
  131. IMPLON=0
  132. SEGACT MCOORD
  133. C LECTURE DU PREMIER POINT . SI C'EST UN ELEMENT ON EN EXTRAIT LE
  134. C DERNIER POINT .
  135. IOB1=0
  136. CALL LIROBJ('MAILLAGE',IOB1,0,IRETOU)
  137. IF (IRETOU.EQ.0) GOTO 1
  138. IPT6=IOB1
  139. SEGACT IPT6
  140. CALL EXTRPO(IOB1,1,IPOIN1)
  141. GOTO 2
  142. 1 CALL LIROBJ('POINT ',IPOIN1,1,IRETOU)
  143. 2 IF (IERR.NE.0) RETURN
  144. C PREMIERE SURFACE
  145. CALL LIRMOT(MCLE,8,ICLE,1)
  146. IF (IERR.NE.0) RETURN
  147. C Y A T-IL DES DENSITES IMPOSEES
  148. 82 IF (IMPOI.EQ.1.OR.ICLE.NE.6) GOTO 80
  149. IMPOI=1
  150. CALL LIRREE(TPOIN1,1,IRETOU)
  151. CALL LIRMOT(MCLE,8,ICLE,1)
  152. IF (IERR.NE.0) RETURN
  153. 80 IF (IMPOF.EQ.1.OR.ICLE.NE.7) GOTO 81
  154. IMPOF=1
  155. CALL LIRREE(TPOIN2,1,IRETOU)
  156. CALL LIRMOT(MCLE,8,ICLE,1)
  157. IF (IERR.NE.0) RETURN
  158. IF (IMPOI.EQ.0) GOTO 82
  159. 81 IF (IMPLON.EQ.1.OR.ICLE.NE.8) GOTO 800
  160. IMPLON=1
  161. CALL LIRREE(DLONG1,1,IRETOU)
  162. CALL LIRMOT(MCLE,8,ICLE,1)
  163. IF (IERR.NE.0) RETURN
  164. goto 82
  165. 800 continue
  166. C ACQUISITION DES DONNEES DE LA PREMIERE SURFACE
  167. GOTO (20,21,22,23,24),ICLE
  168. 20 CALL LIROBJ('POINT ',IPL,1,IRETOU)
  169. IF (IERR.NE.0) RETURN
  170. IREF=4*(IPL-1)
  171. XIPL=XCOOR(IREF+1)
  172. YIPL=XCOOR(IREF+2)
  173. ZIPL=XCOOR(IREF+3)
  174. GOTO 25
  175. 21 CALL LIROBJ('POINT ',IPL,1,IRETOU)
  176. IF (IERR.NE.0) RETURN
  177. IREF=4*(IPL-1)
  178. XIPL=XCOOR(IREF+1)
  179. YIPL=XCOOR(IREF+2)
  180. ZIPL=XCOOR(IREF+3)
  181. GOTO 25
  182. 22 CALL LIROBJ('POINT ',IPL,1,IRETOU)
  183. IF (IERR.NE.0) RETURN
  184. IREF=4*(IPL-1)
  185. XIPL=XCOOR(IREF+1)
  186. YIPL=XCOOR(IREF+2)
  187. ZIPL=XCOOR(IREF+3)
  188. CALL LIROBJ('POINT ',IPL2,1,IRETOU)
  189. IF (IERR.NE.0) RETURN
  190. IREF=4*(IPL2-1)
  191. XIV=XCOOR(IREF+1)-XIPL
  192. YIV=XCOOR(IREF+2)-YIPL
  193. ZIV=XCOOR(IREF+3)-ZIPL
  194. S=SQRT(XIV*XIV+YIV*YIV+ZIV*ZIV)
  195. IF (S.EQ.0.) THEN
  196. CALL ERREUR(21)
  197. RETURN
  198. ENDIF
  199. XIV=XIV/S
  200. YIV=YIV/S
  201. ZIV=ZIV/S
  202. GOTO 25
  203. 23 CALL LIROBJ('POINT ',IPL,1,IRETOU)
  204. IF (IERR.NE.0) RETURN
  205. IREF=4*(IPL-1)
  206. XIPL=XCOOR(IREF+1)
  207. YIPL=XCOOR(IREF+2)
  208. ZIPL=XCOOR(IREF+3)
  209. CALL LIROBJ('POINT ',IPL2,1,IRETOU)
  210. IF (IERR.NE.0) RETURN
  211. IREF=4*(IPL2-1)
  212. XIV=XCOOR(IREF+1)-XIPL
  213. YIV=XCOOR(IREF+2)-YIPL
  214. ZIV=XCOOR(IREF+3)-ZIPL
  215. S=SQRT(XIV**2+YIV**2+ZIV**2)
  216. IF (S.EQ.0.) THEN
  217. CALL ERREUR(21)
  218. RETURN
  219. ENDIF
  220. XIV=XIV/S
  221. YIV=YIV/S
  222. ZIV=ZIV/S
  223. GOTO 25
  224. 24 CALL LIROBJ('POINT ',IPL,1,IRETOU)
  225. IF (IERR.NE.0) RETURN
  226. IREF=4*(IPL-1)
  227. XIPL=XCOOR(IREF+1)
  228. YIPL=XCOOR(IREF+2)
  229. ZIPL=XCOOR(IREF+3)
  230. CALL LIROBJ('POINT ',IAXE,1,IRETOU)
  231. IF (IERR.NE.0) RETURN
  232. IREF=4*(IAXE-1)
  233. XIAXE=XCOOR(IREF+1)-XIPL
  234. YIAXE=XCOOR(IREF+2)-YIPL
  235. ZIAXE=XCOOR(IREF+3)-ZIPL
  236. SIAXE=SQRT(XIAXE**2+YIAXE**2+ZIAXE**2)
  237. IF (SIAXE.EQ.0.) THEN
  238. CALL ERREUR(21)
  239. RETURN
  240. ENDIF
  241. XIAXE=XIAXE/SIAXE
  242. YIAXE=YIAXE/SIAXE
  243. ZIAXE=ZIAXE/SIAXE
  244. CALL LIROBJ('POINT ',IP1,1,IRETOU)
  245. IF (IERR.NE.0) RETURN
  246. IREF=4*(IP1-1)
  247. XIP1=XCOOR(IREF+1)
  248. YIP1=XCOOR(IREF+2)
  249. ZIP1=XCOOR(IREF+3)
  250. GRAYI=SQRT((XIPL-XIP1)**2+(YIPL-YIP1)**2+(ZIPL-ZIP1)**2)
  251. IF (GRAYI.EQ.0.) THEN
  252. CALL ERREUR(21)
  253. RETURN
  254. ENDIF
  255. 25 CONTINUE
  256. C DEUXIEME SURFACE
  257. CALL LIRMOT(MCLE,8,JCLE,1)
  258. IF (IERR.NE.0) RETURN
  259. C Y A T-IL DES DENSITES IMPOSEES
  260. 85 IF (IMPOI.EQ.1.OR.JCLE.NE.6) GOTO 83
  261. IMPOI=1
  262. CALL LIRREE(TPOIN1,1,IRETOU)
  263. CALL LIRMOT(MCLE,8,JCLE,1)
  264. IF (IERR.NE.0) RETURN
  265. 83 IF (IMPOF.EQ.1.OR.JCLE.NE.7) GOTO 84
  266. IMPOF=1
  267. CALL LIRREE(TPOIN2,1,IRETOU)
  268. CALL LIRMOT(MCLE,8,JCLE,1)
  269. IF (IERR.NE.0) RETURN
  270. IF (IMPOI.EQ.0) GOTO 85
  271. 84 IF (IMPLON.EQ.1.OR.ICLE.NE.8) GOTO 801
  272. IMPLON=1
  273. CALL LIRREE(DLONG1,1,IRETOU)
  274. CALL LIRMOT(MCLE,8,ICLE,1)
  275. IF (IERR.NE.0) RETURN
  276. goto 85
  277. 801 continue
  278. C AQUISITION DES DONNEES DE LA DEUXIEME SURFACE
  279. GOTO (40,41,42,43,44),JCLE
  280. 40 CALL LIROBJ('POINT ',JPL,1,IRETOU)
  281. IF (IERR.NE.0) RETURN
  282. IREF=4*(JPL-1)
  283. XJPL=XCOOR(IREF+1)
  284. YJPL=XCOOR(IREF+2)
  285. ZJPL=XCOOR(IREF+3)
  286. GOTO 45
  287. 41 CALL LIROBJ('POINT ',JPL,1,IRETOU)
  288. IF (IERR.NE.0) RETURN
  289. IREF=4*(JPL-1)
  290. XJPL=XCOOR(IREF+1)
  291. YJPL=XCOOR(IREF+2)
  292. ZJPL=XCOOR(IREF+3)
  293. GOTO 45
  294. 42 CALL LIROBJ('POINT ',JPL,1,IRETOU)
  295. IF (IERR.NE.0) RETURN
  296. IREF=4*(JPL-1)
  297. XJPL=XCOOR(IREF+1)
  298. YJPL=XCOOR(IREF+2)
  299. ZJPL=XCOOR(IREF+3)
  300. CALL LIROBJ('POINT ',JPL2,1,IRETOU)
  301. IF (IERR.NE.0) RETURN
  302. IREF=4*(JPL2-1)
  303. XJV=XCOOR(IREF+1)-XJPL
  304. YJV=XCOOR(IREF+2)-YJPL
  305. ZJV=XCOOR(IREF+3)-ZJPL
  306. S=SQRT(XJV**2+YJV**2+ZJV**2)
  307. IF(S.EQ.0.) THEN
  308. CALL ERREUR(21)
  309. RETURN
  310. ENDIF
  311. XJV=XJV/S
  312. YJV=YJV/S
  313. ZJV=ZJV/S
  314. GOTO 45
  315. 43 CALL LIROBJ('POINT ',JPL,1,IRETOU)
  316. IF (IERR.NE.0) RETURN
  317. IREF=4*(JPL-1)
  318. XJPL=XCOOR(IREF+1)
  319. YJPL=XCOOR(IREF+2)
  320. ZJPL=XCOOR(IREF+3)
  321. CALL LIROBJ('POINT ',JPL2,1,IRETOU)
  322. IF (IERR.NE.0) RETURN
  323. IREF=4*(JPL2-1)
  324. XJV=XCOOR(IREF+1)-XJPL
  325. YJV=XCOOR(IREF+2)-YJPL
  326. ZJV=XCOOR(IREF+3)-ZJPL
  327. S=SQRT(XJV**2+YJV**2+ZJV**2)
  328. IF (S.EQ.0.) THEN
  329. CALL ERREUR(21)
  330. RETURN
  331. ENDIF
  332. XJV=XJV/S
  333. YJV=YJV/S
  334. ZJV=ZJV/S
  335. GOTO 45
  336. 44 CALL LIROBJ('POINT ',JPL,1,IRETOU)
  337. IF (IERR.NE.0) RETURN
  338. IREF=4*(JPL-1)
  339. XJPL=XCOOR(IREF+1)
  340. YJPL=XCOOR(IREF+2)
  341. ZJPL=XCOOR(IREF+3)
  342. CALL LIROBJ('POINT ',JAXE,1,IRETOU)
  343. IF (IERR.NE.0) RETURN
  344. IREF=4*(JAXE-1)
  345. XJAXE=XCOOR(IREF+1)-XJPL
  346. YJAXE=XCOOR(IREF+2)-YJPL
  347. ZJAXE=XCOOR(IREF+3)-ZJPL
  348. SJAXE=SQRT(XJAXE**2+YJAXE**2+ZJAXE**2)
  349. IF (SJAXE.EQ.0.) THEN
  350. CALL ERREUR(21)
  351. RETURN
  352. ENDIF
  353. XJAXE=XJAXE/SJAXE
  354. YJAXE=YJAXE/SJAXE
  355. ZJAXE=ZJAXE/SJAXE
  356. CALL LIROBJ('POINT ',JP1,1,IRETOU)
  357. IF (IERR.NE.0) RETURN
  358. IREF=4*(JP1-1)
  359. XJP1=XCOOR(IREF+1)
  360. YJP1=XCOOR(IREF+2)
  361. ZJP1=XCOOR(IREF+3)
  362. GRAYJ=SQRT((XJPL-XJP1)**2+(YJPL-YJP1)**2+(ZJPL-ZJP1)**2)
  363. IF (GRAYJ.EQ.0.) THEN
  364. CALL ERREUR(21)
  365. RETURN
  366. ENDIF
  367. 45 CONTINUE
  368. C Y A T-IL DES DENSITES IMPOSEES
  369. IRETOU=0
  370. 88 IF (IMPOI*IMPOF*implon.EQ.0) CALL LIRMOT(MCLE(6),3,MDPOS,0)
  371. IF (MDPOS.EQ.0) GOTO 89
  372. IF (IMPOI.EQ.1.OR.MDPOS.NE.1) GOTO 86
  373. IMPOI=1
  374. CALL LIRREE(TPOIN1,1,IRETOU)
  375. goto 88
  376. 86 IF (IMPOF.EQ.1.OR.MDPOS.NE.2) GOTO 87
  377. IMPOF=1
  378. CALL LIRREE(TPOIN2,1,IRETOU)
  379. GOTO 88
  380. 87 if (implon.EQ.1.OR.MDPOS.NE.3) GOTO 802
  381. implon=1
  382. CALL LIRREE(dlong1,1,IRETOU)
  383. goto 88
  384. 802 CALL REFUS
  385. 89 CONTINUE
  386. if (ierr.ne.0) return
  387. C LECTURE DU DEUXIEME POINT . SI C'EST UN ELEMENT ON EN EXTRAIT LE
  388. C PREMIER POINT .
  389. IOB2=0
  390. iretou=0
  391. if (implon.eq.0) CALL LIROBJ('MAILLAGE',IOB2,0,IRETOU)
  392. IF (IRETOU.EQ.0) GOTO 50
  393. IPT6=IOB2
  394. SEGACT IPT6
  395. CALL EXTRPO(IOB2,2,IPOIN2)
  396. GOTO 51
  397. 50 CALL LIROBJ('POINT ',IPOIN2,1,IRETOU)
  398. 51 IF (IERR.NE.0) RETURN
  399. NBP=NBT-2
  400. IREF1=(IPOIN1-1)*4
  401. IREF2=(IPOIN2-1)*4
  402. XPOIN1=XCOOR(IREF1+1)
  403. YPOIN1=XCOOR(IREF1+2)
  404. ZPOIN1=XCOOR(IREF1+3)
  405. IF (IMPOI.EQ.0) TPOIN1=XCOOR(IREF1+4)
  406. XPOIN2=XCOOR(IREF2+1)
  407. YPOIN2=XCOOR(IREF2+2)
  408. ZPOIN2=XCOOR(IREF2+3)
  409. IF (IMPOF.EQ.0) TPOIN2=XCOOR(IREF2+4)
  410. XVEC=XPOIN2-XPOIN1
  411. YVEC=YPOIN2-YPOIN1
  412. ZVEC=ZPOIN2-ZPOIN1
  413. SVEC=SQRT(XVEC**2+YVEC**2+ZVEC**2)
  414. IF (SVEC.EQ.0.) THEN
  415. CALL ERREUR(21)
  416. RETURN
  417. ENDIF
  418. EPS2=1E-10*SVEC**2
  419. C IL FAUT VERIFIER LA COMPATIBILITE DES DONNEES
  420. GOTO (60,61,62,63,64),ICLE
  421. 60 CONTINUE
  422. XPV1=XPOIN1-XIPL
  423. YPV1=YPOIN1-YIPL
  424. ZPV1=ZPOIN1-ZIPL
  425. XPV2=XPOIN2-XIPL
  426. YPV2=YPOIN2-YIPL
  427. ZPV2=ZPOIN2-ZIPL
  428. XIPN=YPV1*ZPV2-ZPV1*YPV2
  429. YIPN=ZPV1*XPV2-XPV1*ZPV2
  430. ZIPN=XPV1*YPV2-YPV1*XPV2
  431. SIPN=SQRT(XIPN**2+YIPN**2+ZIPN**2)
  432. IF (SIPN.EQ.0.) CALL ERREUR(21)
  433. IF (IERR.NE.0) RETURN
  434. XIPN=XIPN/SIPN
  435. YIPN=YIPN/SIPN
  436. ZIPN=ZIPN/SIPN
  437. GOTO 65
  438. 61 CONTINUE
  439. RAYI=SQRT((XPOIN1-XIPL)**2+(YPOIN1-YIPL)**2+(ZPOIN1-ZIPL)**2)
  440. RAYB=SQRT((XPOIN2-XIPL)**2+(YPOIN2-YIPL)**2+(ZPOIN2-ZIPL)**2)
  441. if (implon.eq.0) then
  442. IF (ABS((RAYI-RAYB)/RAYI).GE.1E-2) THEN
  443. CALL ERREUR(21)
  444. RETURN
  445. ENDIF
  446. RAYI=SQRT(RAYI*RAYB)
  447. endif
  448. GOTO 65
  449. 62 CONTINUE
  450. U=XPOIN1-XIPL
  451. V=YPOIN1-YIPL
  452. W=ZPOIN1-ZIPL
  453. SCA=U*XIV+V*YIV+W*ZIV
  454. U=U-SCA*XIV
  455. V=V-SCA*YIV
  456. W=W-SCA*ZIV
  457. RAYI=SQRT(U**2+V**2+W**2)
  458. U=XPOIN2-XIPL
  459. V=YPOIN2-YIPL
  460. W=ZPOIN2-ZIPL
  461. SCA=U*XIV+V*YIV+W*ZIV
  462. U=U-SCA*XIV
  463. V=V-SCA*YIV
  464. W=W-SCA*ZIV
  465. RAYB=SQRT(U**2+V**2+W**2)
  466. if (implon.eq.0) then
  467. IF (ABS(RAYI-RAYB)/RAYI.GT.1E-2) CALL ERREUR(21)
  468. RAYI=SQRT(RAYI*RAYB)
  469. endif
  470. GOTO 65
  471. 63 CONTINUE
  472. U=XPOIN1-XIPL
  473. V=YPOIN1-YIPL
  474. W=ZPOIN1-ZIPL
  475. SCA=ABS(U*XIV+V*YIV+W*ZIV)
  476. U1=V*ZIV-W*YIV
  477. V1=W*XIV-U*ZIV
  478. W1=U*YIV-V*XIV
  479. VEC=SQRT(U1**2+V1**2+W1**2)
  480. IF (SCA.EQ.0.) CALL ERREUR(21)
  481. IF (IERR.NE.0) RETURN
  482. TANGI=VEC/SCA
  483. U=XPOIN2-XIPL
  484. V=YPOIN2-YIPL
  485. W=ZPOIN2-ZIPL
  486. SCA=ABS(U*XIV+V*YIV+W*ZIV)
  487. U1=V*ZIV-W*YIV
  488. V1=W*XIV-U*ZIV
  489. W1=U*YIV-V*XIV
  490. VEC=SQRT(U1**2+V1**2+W1**2)
  491. IF (SCA.EQ.0.) CALL ERREUR(21)
  492. IF (IERR.NE.0) RETURN
  493. TANGB=VEC/SCA
  494. if (implon.eq.0) then
  495. IF (ABS(TANGI-TANGB)/TANGI.GE.1E-2) CALL ERREUR(21)
  496. TANGI=SQRT(TANGI*TANGB)
  497. endif
  498. GOTO 65
  499. 64 CONTINUE
  500. U=XPOIN1-XIPL
  501. V=YPOIN1-YIPL
  502. W=ZPOIN1-ZIPL
  503. SCA=U*XIAXE+V*YIAXE+W*ZIAXE
  504. U=U-SCA*XIAXE
  505. V=V-SCA*YIAXE
  506. W=W-SCA*ZIAXE
  507. S=SQRT(U**2+V**2+W**2)
  508. IF (S.EQ.0.) CALL ERREUR(21)
  509. IF (IERR.NE.0) RETURN
  510. XC=XIPL+U*GRAYI/S
  511. YC=YIPL+V*GRAYI/S
  512. ZC=ZIPL+W*GRAYI/S
  513. PRAYI=SQRT((XPOIN1-XC)**2+(YPOIN1-YC)**2+(ZPOIN1-ZC)**2)
  514. U=XPOIN2-XIPL
  515. V=YPOIN2-YIPL
  516. W=ZPOIN2-ZIPL
  517. SCA=U*XIAXE+V*YIAXE+W*ZIAXE
  518. U=U-SCA*XIAXE
  519. V=V-SCA*YIAXE
  520. W=W-SCA*ZIAXE
  521. S=SQRT(U**2+V**2+W**2)
  522. IF (S.EQ.0.) CALL ERREUR(21)
  523. IF (IERR.NE.0) RETURN
  524. XC=XIPL+U*GRAYI/S
  525. YC=YIPL+V*GRAYI/S
  526. ZC=ZIPL+W*GRAYI/S
  527. PRAYB=SQRT((XPOIN2-XC)**2+(YPOIN2-YC)**2+(ZPOIN2-ZC)**2)
  528. if (implon.eq.0) then
  529. IF (ABS(PRAYI-PRAYB)/PRAYI.GE.1E-2) CALL ERREUR(21)
  530. PRAYI=SQRT(PRAYI*PRAYB)
  531. endif
  532. 65 IF (IERR.NE.0) RETURN
  533. GOTO (70,71,72,73,74),JCLE
  534. 70 CONTINUE
  535. XPV1=XPOIN1-XJPL
  536. YPV1=YPOIN1-YJPL
  537. ZPV1=ZPOIN1-ZJPL
  538. XPV2=XPOIN2-XJPL
  539. YPV2=YPOIN2-YJPL
  540. ZPV2=ZPOIN2-ZJPL
  541. XJPN=YPV1*ZPV2-ZPV1*YPV2
  542. YJPN=ZPV1*XPV2-XPV1*ZPV2
  543. ZJPN=XPV1*YPV2-YPV1*XPV2
  544. SJPN=SQRT(XJPN**2+YJPN**2+ZJPN**2)
  545. IF (SJPN.EQ.0.) CALL ERREUR(21)
  546. IF (IERR.NE.0) RETURN
  547. XJPN=XJPN/SJPN
  548. YJPN=YJPN/SJPN
  549. ZJPN=ZJPN/SJPN
  550. GOTO 75
  551. 71 CONTINUE
  552. RAYJ=SQRT((XPOIN1-XJPL)**2+(YPOIN1-YJPL)**2+(ZPOIN1-ZJPL)**2)
  553. RAYB=SQRT((XPOIN2-XJPL)**2+(YPOIN2-YJPL)**2+(ZPOIN2-ZJPL)**2)
  554. if (implon.eq.0) then
  555. IF (ABS((RAYJ-RAYB)/RAYJ).GE.1E-2) CALL ERREUR(21)
  556. RAYJ=SQRT(RAYJ*RAYB)
  557. endif
  558. GOTO 75
  559. 72 CONTINUE
  560. U=XPOIN1-XJPL
  561. V=YPOIN1-YJPL
  562. W=ZPOIN1-ZJPL
  563. SCA=U*XJV+V*YJV+W*ZJV
  564. U=U-SCA*XJV
  565. V=V-SCA*YJV
  566. W=W-SCA*ZJV
  567. RAYJ=SQRT(U**2+V**2+W**2)
  568. U=XPOIN2-XJPL
  569. V=YPOIN2-YJPL
  570. W=ZPOIN2-ZJPL
  571. SCA=U*XJV+V*YJV+W*ZJV
  572. U=U-SCA*XJV
  573. V=V-SCA*YJV
  574. W=W-SCA*ZJV
  575. RAYB=SQRT(U**2+V**2+W**2)
  576. if (implon.eq.0) then
  577. IF (ABS(RAYJ-RAYB)/RAYJ.GT.1E-2) CALL ERREUR(21)
  578. RAYJ=SQRT(RAYJ*RAYB)
  579. endif
  580. GOTO 75
  581. 73 CONTINUE
  582. U=XPOIN1-XJPL
  583. V=YPOIN1-YJPL
  584. W=ZPOIN1-ZJPL
  585. SCA=ABS(U*XJV+V*YJV+W*ZJV)
  586. U1=V*ZJV-W*YJV
  587. V1=W*XJV-U*ZJV
  588. W1=U*YJV-V*XJV
  589. VEC=SQRT(U1**2+V1**2+W1**2)
  590. IF (SCA.EQ.0.) CALL ERREUR(21)
  591. IF (IERR.NE.0) RETURN
  592. TANGJ=VEC/SCA
  593. U=XPOIN2-XJPL
  594. V=YPOIN2-YJPL
  595. W=ZPOIN2-ZJPL
  596. SCA=ABS(U*XJV+V*YJV+W*ZJV)
  597. U1=V*ZJV-W*YJV
  598. V1=W*XJV-U*ZJV
  599. W1=U*YJV-V*XJV
  600. VEC=SQRT(U1**2+V1**2+W1**2)
  601. IF (SCA.EQ.0.) CALL ERREUR(21)
  602. IF (IERR.NE.0) RETURN
  603. TANGB=VEC/SCA
  604. if (implon.eq.0) then
  605. IF (ABS(TANGJ-TANGB)/TANGJ.GE.1E-2) CALL ERREUR(21)
  606. TANGJ=SQRT(TANGJ*TANGB)
  607. endif
  608. GOTO 75
  609. 74 CONTINUE
  610. U=XPOIN1-XJPL
  611. V=YPOIN1-YJPL
  612. W=ZPOIN1-ZJPL
  613. SCA=U*XJAXE+V*YJAXE+W*ZJAXE
  614. U=U-SCA*XJAXE
  615. V=V-SCA*YJAXE
  616. W=W-SCA*ZJAXE
  617. S=SQRT(U**2+V**2+W**2)
  618. IF (S.EQ.0.) CALL ERREUR(21)
  619. IF (IERR.NE.0) RETURN
  620. XC=XJPL+U*GRAYJ/S
  621. YC=YJPL+V*GRAYJ/S
  622. ZC=ZJPL+W*GRAYJ/S
  623. PRAYJ=SQRT((XPOIN1-XC)**2+(YPOIN1-YC)**2+(ZPOIN1-ZC)**2)
  624. U=XPOIN2-XJPL
  625. V=YPOIN2-YJPL
  626. W=ZPOIN2-ZJPL
  627. SCA=U*XJAXE+V*YJAXE+W*ZJAXE
  628. U=U-SCA*XJAXE
  629. V=V-SCA*YJAXE
  630. W=W-SCA*ZJAXE
  631. S=SQRT(U**2+V**2+W**2)
  632. IF (S.EQ.0.) CALL ERREUR(21)
  633. IF (IERR.NE.0) RETURN
  634. XC=XJPL+U*GRAYJ/S
  635. YC=YJPL+V*GRAYJ/S
  636. ZC=ZJPL+W*GRAYJ/S
  637. PRAYB=SQRT((XPOIN2-XC)**2+(YPOIN2-YC)**2+(ZPOIN2-ZC)**2)
  638. if (implon.eq.0) then
  639. IF (ABS(PRAYJ-PRAYB)/PRAYJ.GE.1E-2) CALL ERREUR(21)
  640. PRAYJ=SQRT(PRAYJ*PRAYB)
  641. endif
  642. 75 IF (IERR.NE.0) RETURN
  643. C ON PLACE NBP POINTS SUR L'INTERSECTION POUR CALCULER DES ABCISSES
  644. C CURVILIGNES
  645. XVECR=XVEC/(NBP+1)
  646. YVECR=YVEC/(NBP+1)
  647. ZVECR=ZVEC/(NBP+1)
  648. XVEC=XVEC/SVEC
  649. YVEC=YVEC/SVEC
  650. ZVEC=ZVEC/SVEC
  651. XPAUX(1,1)=XPOIN1
  652. XPAUX(2,1)=YPOIN1
  653. XPAUX(3,1)=ZPOIN1
  654. XPAUX(4,1)=0.
  655. IBP=0
  656. 102 IBP=IBP+1
  657. IF (IBP.GT.NBP) GOTO 100
  658. XP=XPAUX(1,IBP)+XVECR
  659. YP=XPAUX(2,IBP)+YVECR
  660. ZP=XPAUX(3,IBP)+ZVECR
  661. XPA=XP
  662. YPA=YP
  663. ZPA=ZP
  664. IRET = 101
  665. GOTO 400
  666. 101 CONTINUE
  667. IF (IIMPI.EQ.1) WRITE (IOIMP,8880) IBP,XPA,YPA,ZPA,XP1,YP1,ZP1
  668. 8880 FORMAT(I4,' POINT A PROJETER ',3G13.5,' POINT PROJETE ',3G13.5)
  669. C ON RANGE LE POINT D'INTERSECTION
  670. XPAUX(1,IBP+1)=XP1
  671. XPAUX(2,IBP+1)=YP1
  672. XPAUX(3,IBP+1)=ZP1
  673. XPAUX(4,IBP+1)=XPAUX(4,IBP)+SQRT((XP1-XPAUX(1,IBP))**2+
  674. # (YP1-XPAUX(2,IBP))**2+(ZP1-XPAUX(3,IBP))**2)
  675. GOTO 102
  676. 100 CONTINUE
  677. XPAUX(1,NBT)=XPOIN2
  678. XPAUX(2,NBT)=YPOIN2
  679. XPAUX(3,NBT)=ZPOIN2
  680. XPAUX(4,NBT)=XPAUX(4,NBT-1)+SQRT((XP1-XPOIN2)**2+(YP1-YPOIN2)**2+
  681. # (ZP1-ZPOIN2)**2)
  682. DLONG=XPAUX(4,NBT)
  683. if (implon.eq.1) then
  684. if ( iimpi.eq.1) write(ioimp,8888) dlong1 ,dlong
  685. 8888 format ( ' longueur demandee et longueur donnee',2g13.5)
  686. if (dlong1.gt.dlong) call erreur(21)
  687. if (ierr.ne.0) return
  688. dlong=dlong1
  689. endif
  690. DEN1=TPOIN1/DLONG
  691. DEN2=TPOIN2/DLONG
  692. CALL DECOUP(INBR,DEN1,DEN2,APROG,NBELEM,DENI,DECA,DLONG)
  693. NX=NBELEM-1
  694. DEN1=DEN1*DLONG
  695. IF (IIMPI.EQ.1) WRITE (IOIMP,9900) DLONG,NBELEM,APROG
  696. 9900 FORMAT (' LONGUEUR ',G12.5,' ELEMENT ',I5,' RAISON ',G12.5)
  697. C CREATION DU SEGMENT
  698. NBNN=NBNNE(KDEGRE(ILCOUR))
  699. IF (ILCOUR.EQ.0) CALL ERREUR(16)
  700. IF (KDEGRE(ILCOUR).EQ.0) CALL ERREUR(16)
  701. IF (NBNN.EQ.0) CALL ERREUR(16)
  702. IF (IERR.NE.0) RETURN
  703. NBSOUS=0
  704. NBREF=0
  705. SEGINI MELEME
  706. ITYPEL=KDEGRE(ILCOUR)
  707. NPOIN=(ITYPEL-1)*NBELEM-1
  708. if (implon.eq.1) npoin=npoin+1
  709. SEGINI TABPAR
  710. C DEFINITION DE LA TOPOLOGIE
  711. IPOO=XCOOR(/1)/(IDIM+1)
  712. IADR=IPOO
  713. DIN=DEN1
  714. DL=0.
  715. NUM(1,1)=IPOIN1
  716. IF (NX.EQ.0) GOTO 510
  717. DO 512 I=1,NX
  718. DIN=DIN*APROG
  719. IF (ITYPEL.EQ.2) GOTO 513
  720. IPOO=IPOO+1
  721. NUM(2,I)=IPOO
  722. TABPAR(IPOO-IADR)=DIN/2.+DL
  723. 513 IPOO=IPOO+1
  724. NUM(ITYPEL,I)=IPOO
  725. NUM(1,I+1)=IPOO
  726. DL=DL+DIN
  727. TABPAR(IPOO-IADR)=DL
  728. 512 CONTINUE
  729. 510 CONTINUE
  730. IF (ITYPEL.NE.3) GOTO 514
  731. DIN=DIN*APROG
  732. IPOO=IPOO+1
  733. TABPAR(IPOO-IADR)=DL+DIN/2
  734. NUM(2,NBELEM)=IPOO
  735. 514 NUM(ITYPEL,NBELEM)=IPOIN2
  736. if (implon.eq.1) then
  737. NUM(ITYPEL,NBELEM)=ipoo+1
  738. TABPAR(IPOO+1-IADR)=DLong
  739. endif
  740. C CREATION DES POINTS LA LONGUEUR A RESPECTER EST DANS TABPAR
  741. C LA LONGUEUR DES POINTS CONNUS EST DANS XPAUX
  742. IP=1
  743. T=0.
  744. DLR=0.
  745. NBPTA=XCOOR(/1)/(IDIM+1)
  746. NBPTS=NBPTA+NPOIN
  747. SEGADJ MCOORD
  748. I=0
  749. 523 I=I+1
  750. IF (I.GT.NPOIN) GOTO 520
  751. DLV=TABPAR(I)
  752. 522 IF (DLV.LE.DLR) GOTO 521
  753. IP=IP+1
  754. DLR=XPAUX(4,IP)
  755. GOTO 522
  756. 521 CONTINUE
  757. DLA=XPAUX(4,IP-1)
  758. COF=(DLV-DLA)/(DLR-DLA)
  759. XP=XPAUX(1,IP-1)+COF*(XPAUX(1,IP)-XPAUX(1,IP-1))
  760. YP=XPAUX(2,IP-1)+COF*(XPAUX(2,IP)-XPAUX(2,IP-1))
  761. ZP=XPAUX(3,IP-1)+COF*(XPAUX(3,IP)-XPAUX(3,IP-1))
  762. IF (IIMPI.NE.0) WRITE (IOIMP,8878) I,IP,XP,YP,ZP
  763. 8878 FORMAT(' AVANT ',I4,' INTERPOLE ',I4,3G13.5)
  764. IRET = 525
  765. GOTO 400
  766. 525 CONTINUE
  767. IF (IIMPI.NE.0) WRITE (IOIMP,8876) I,IP,XP1,YP1,ZP1
  768. 8876 FORMAT(' APRES ',I4,' INTERPOLE ',I4,3G13.5)
  769. XCOOR(NBPTA*4+1)=XP1
  770. XCOOR(NBPTA*4+2)=YP1
  771. XCOOR(NBPTA*4+3)=ZP1
  772. XCOOR(NBPTA*4+4)=TABPAR(I)-T
  773. NBPTA=NBPTA+1
  774. T=TABPAR(I)
  775. IADR=IADR+1
  776. GOTO 523
  777. 520 CONTINUE
  778. SEGSUP TABPAR
  779. IPT2=MELEME
  780. DO 527 I=1,NUM(/2)
  781. 527 ICOLOR(I)=IDCOUL
  782. IF (IOB1.NE.0) THEN
  783. ltelq=.false.
  784. CALL FUSE(IOB1,MELEME,IPT2,ltelq)
  785. IPT6=IOB1
  786. SEGDES IPT6
  787. ENDIF
  788. MELEME=IPT2
  789. IF (IOB2.NE.0) THEN
  790. ltelq=.false.
  791. CALL FUSE(MELEME,IOB2,IPT2,ltelq)
  792. IPT6=IOB2
  793. SEGDES IPT6
  794. ENDIF
  795. MELEME=IPT2
  796. SEGDES MELEME
  797. CALL ECROBJ('MAILLAGE',MELEME)
  798. RETURN
  799. C SOUS-PROGRAMME DE CALCUL ITERATIF D'INTERSECTION
  800. 400 CONTINUE
  801. ICYCL=0
  802. 402 CONTINUE
  803. C IL FAUT TROUVER L'INTERSECTION DU PLAN PERPENDICULAIRE A XVEC YVEC
  804. C ZVEC EN XPOIN YPOIN ZPOIN ET DES DEUX SURFACES
  805. C ON UTILISE UNE METHODE ITERATIVE
  806. C ON VA CHERCHER LE POINT DE PROJECTION DE XP SUR CHAQUE SURFACE DANS
  807. C LE PLAN ET UN VECTEUR NORMAL
  808. C XP YP ZP POINT A PROJETER
  809. C XP1 YP1 ZP1 INTERSECTION EN RETOUR
  810. GOTO (200,210,220,230,240),ICLE
  811. 250 CONTINUE
  812. GOTO (300,310,320,330,340),JCLE
  813. 350 CONTINUE
  814. C LES PTS ET LES VECTS SONT XPI... XPJ... XNI... XNJ...
  815. C ON CALCULE L'INTERSECTION DES DEUX VECTEURS (AFFINES)
  816. XDIF=XPJ-XPI
  817. YDIF=YPJ-YPI
  818. ZDIF=ZPJ-ZPI
  819. C COMPOSANTES SUR CHAQUE VECTEUR
  820. DET=(YNI*ZNJ-ZNI*YNJ)*XVEC
  821. # +(ZNI*XNJ-XNI*ZNJ)*YVEC
  822. # +(XNI*YNJ-YNI*XNJ)*ZVEC
  823. IF (DET**2.LT.EPS2) CALL ERREUR(40)
  824. A=((YDIF*ZNJ-ZDIF*YNJ)*XVEC
  825. # +(ZDIF*XNJ-XDIF*ZNJ)*YVEC
  826. # +(XDIF*YNJ-YDIF*XNJ)*ZVEC)/DET
  827. XP1=XPI+A*XNI
  828. YP1=YPI+A*YNI
  829. ZP1=ZPI+A*ZNI
  830. C TEST CONVERGENCE
  831. ECART=(XP-XP1)**2+(YP-YP1)**2+(ZP-ZP1)**2
  832. IF (ECART.LE.EPS2) THEN
  833. IF (IRET.EQ.101) GOTO 101
  834. IF (IRET.EQ.525) GOTO 525
  835. ENDIF
  836. XP=XP1
  837. YP=YP1
  838. ZP=ZP1
  839. ICYCL=ICYCL+1
  840. IF (IIMPI.EQ.1) WRITE (IOIMP,4478) ICYCL,ECART,EPS2,XPI,YPI,ZPI,
  841. # XNI,YNI,ZNI,XPJ,YPJ,ZPJ,XNJ,YNJ,ZNJ,XP,YP,ZP
  842. 4478 FORMAT(' ICYCL ',I4,' ECART ',G12.5,' EPS2 ',G12.5,/,' PI ',
  843. # 3G13.5,' NI ',3G13.5,/,' PJ ',3G13.5,' NJ ',3G13.5,' P ',3G13.5)
  844. IF (ICYCL.LE.25) GOTO 402
  845. CALL ERREUR(40)
  846. RETURN
  847. C PROJECTION SUR UN PLAN XIPL XIPN EN RESTANT DANS XVEC
  848. 200 CONTINUE
  849. C VECT PERPENDICULAIRE A XVEC ET XIPN
  850. U=YVEC*ZIPN-ZVEC*YIPN
  851. V=ZVEC*XIPN-XVEC*ZIPN
  852. W=XVEC*YIPN-YVEC*XIPN
  853. S=SQRT(U**2+V**2+W**2)
  854. IF (S.EQ.0.) CALL ERREUR(40)
  855. U1=U/S
  856. V1=V/S
  857. W1=W/S
  858. U=V1*ZVEC-W1*YVEC
  859. V=W1*XVEC-U1*ZVEC
  860. W=U1*YVEC-V1*XVEC
  861. C VECT SUIVANT LEQUEL ON DEPLACE LE POINT
  862. PAR=-((XP-XIPL)*XIPN+(YP-YIPL)*YIPN+(ZP-ZIPL)*ZIPN)/
  863. # (U*XIPN+V*YIPN+W*ZIPN)
  864. XPI=XP+PAR*U
  865. YPI=YP+PAR*V
  866. ZPI=ZP+PAR*W
  867. XNI=U1
  868. YNI=V1
  869. ZNI=W1
  870. GOTO 250
  871. 300 CONTINUE
  872. U=YVEC*ZJPN-ZVEC*YJPN
  873. V=ZVEC*XJPN-XVEC*ZJPN
  874. W=XVEC*YJPN-YVEC*XJPN
  875. S=SQRT(U**2+V**2+W**2)
  876. IF (S.EQ.0.) CALL ERREUR(40)
  877. U1=U/S
  878. V1=V/S
  879. W1=W/S
  880. U=V1*ZVEC-W1*YVEC
  881. V=W1*XVEC-U1*ZVEC
  882. W=U1*YVEC-V1*XVEC
  883. PAR=-((XP-XJPL)*XJPN+(YP-YJPL)*YJPN+(ZP-ZJPL)*ZJPN)/
  884. # (U*XJPN+V*YJPN+W*ZJPN)
  885. XPJ=XP+PAR*U
  886. YPJ=YP+PAR*V
  887. ZPJ=ZP+PAR*W
  888. XNJ=U1
  889. YNJ=V1
  890. ZNJ=W1
  891. GOTO 350
  892. C PROJECTION SUR SPHERE CENTRE XIPL RAYON RAYI EN RESTANT SUR LE PLAN X
  893. 210 CONTINUE
  894. U1=XIPL-XP
  895. V1=YIPL-YP
  896. W1=ZIPL-ZP
  897. SCA=U1*XVEC+V1*YVEC+W1*ZVEC
  898. XC=XIPL-SCA*XVEC
  899. YC=YIPL-SCA*YVEC
  900. ZC=ZIPL-SCA*ZVEC
  901. IF (SCA.GT.RAYI) CALL ERREUR(40)
  902. IF (IERR.NE.0) RETURN
  903. RAYC=SQRT(RAYI**2-SCA**2)
  904. U=XC-XP
  905. V=YC-YP
  906. W=ZC-ZP
  907. S=SQRT(U**2+V**2+W**2)
  908. IF (S.EQ.0.) CALL ERREUR(40)
  909. IF (IERR.NE.0) RETURN
  910. U=U/S
  911. V=V/S
  912. W=W/S
  913. XPI=XC-RAYC*U
  914. YPI=YC-RAYC*V
  915. ZPI=ZC-RAYC*W
  916. XNI=YVEC*W-ZVEC*V
  917. YNI=ZVEC*U-XVEC*W
  918. ZNI=XVEC*V-YVEC*U
  919. GOTO 250
  920. 310 CONTINUE
  921. U1=XJPL-XP
  922. V1=YJPL-YP
  923. W1=ZJPL-ZP
  924. SCA=U1*XVEC+V1*YVEC+W1*ZVEC
  925. XC=XJPL-SCA*XVEC
  926. YC=YJPL-SCA*YVEC
  927. ZC=ZJPL-SCA*ZVEC
  928. IF (SCA.GT.RAYJ) CALL ERREUR(40)
  929. IF (IERR.NE.0) RETURN
  930. RAYC=SQRT(RAYJ**2-SCA**2)
  931. U=XC-XP
  932. V=YC-YP
  933. W=ZC-ZP
  934. S=SQRT(U**2+V**2+W**2)
  935. IF (S.EQ.0.) CALL ERREUR(40)
  936. IF (IERR.NE.0) RETURN
  937. U=U/S
  938. V=V/S
  939. W=W/S
  940. XPJ=XC-RAYC*U
  941. YPJ=YC-RAYC*V
  942. ZPJ=ZC-RAYC*W
  943. XNJ=YVEC*W-ZVEC*V
  944. YNJ=ZVEC*U-XVEC*W
  945. ZNJ=XVEC*V-YVEC*U
  946. GOTO 350
  947. C PROJECTION SUR CYLINDRE AXE XIPL DIR XIV EN RESTANT DANS XVEC
  948. 220 CONTINUE
  949. U=XIPL-XP
  950. V=YIPL-YP
  951. W=ZIPL-ZP
  952. SCA=U*XIV+V*YIV+W*ZIV
  953. U=U-SCA*XIV
  954. V=V-SCA*YIV
  955. W=W-SCA*ZIV
  956. S=SQRT(U**2+V**2+W**2)
  957. IF (S.EQ.0.) CALL ERREUR(40)
  958. IF (IERR.NE.0) RETURN
  959. U1=U-U*RAYI/S
  960. V1=V-V*RAYI/S
  961. W1=W-W*RAYI/S
  962. DNUM=U1*XVEC+V1*YVEC+W1*ZVEC
  963. DNOM=XIV*XVEC+YIV*YVEC+ZIV*ZVEC
  964. IF (DNOM.EQ.0.) GOTO 221
  965. IF (IERR.NE.0) RETURN
  966. DPAR=-DNUM/DNOM
  967. XPI=XP+U1+DPAR*XIV
  968. YPI=YP+V1+DPAR*YIV
  969. ZPI=ZP+W1+DPAR*ZIV
  970. XNI=(V*ZVEC-W*YVEC)/S
  971. YNI=(W*XVEC-U*ZVEC)/S
  972. ZNI=(U*YVEC-V*XVEC)/S
  973. GOTO 250
  974. C L'AXE DU CYLINDRE EST PARALLELE AU PLAN XVEC
  975. 221 SCA=U*XVEC+V*YVEC+W*ZVEC
  976. IF (SCA.GT.RAYI) CALL ERREUR(40)
  977. IF (IERR.NE.0) RETURN
  978. RAYC=SQRT(RAYI**2-SCA**2)
  979. U=U-SCA*XVEC
  980. V=V-SCA*YVEC
  981. W=W-SCA*ZVEC
  982. S=SQRT(U**2+V**2+W**2)
  983. IF (S.EQ.0.) CALL ERREUR(40)
  984. IF (IERR.NE.0) RETURN
  985. XPI=XP+U-U*RAYC/S
  986. YPI=YP+V-V*RAYC/S
  987. ZPI=ZP+W-W*RAYC/S
  988. XNI=XIV
  989. YNI=YIV
  990. ZNI=ZIV
  991. GOTO 250
  992. 320 CONTINUE
  993. U=XJPL-XP
  994. V=YJPL-YP
  995. W=ZJPL-ZP
  996. SCA=U*XJV+V*YJV+W*ZJV
  997. U=U-SCA*XJV
  998. V=V-SCA*YJV
  999. W=W-SCA*ZJV
  1000. S=SQRT(U**2+V**2+W**2)
  1001. IF (S.EQ.0.) CALL ERREUR(40)
  1002. IF (IERR.NE.0) RETURN
  1003. U1=U-U*RAYJ/S
  1004. V1=V-V*RAYJ/S
  1005. W1=W-W*RAYJ/S
  1006. DNUM=U1*XVEC+V1*YVEC+W1*ZVEC
  1007. DNOM=XJV*XVEC+YJV*YVEC+ZJV*ZVEC
  1008. IF (DNOM.EQ.0.) GOTO 321
  1009. IF (IERR.NE.0) RETURN
  1010. DPAR=-DNUM/DNOM
  1011. XPJ=XP+U1+DPAR*XJV
  1012. YPJ=YP+V1+DPAR*YJV
  1013. ZPJ=ZP+W1+DPAR*ZJV
  1014. XNJ=(V*ZVEC-W*YVEC)/S
  1015. YNJ=(W*XVEC-U*ZVEC)/S
  1016. ZNJ=(U*YVEC-V*XVEC)/S
  1017. GOTO 350
  1018. C L'AXE DU CYLINDRE EST PARALLELE AU PLAN XVEC
  1019. 321 SCA=U*XVEC+V*YVEC+W*ZVEC
  1020. IF (SCA.GT.RAYJ) CALL ERREUR(40)
  1021. IF (IERR.NE.0) RETURN
  1022. RAYC=SQRT(RAYJ**2-SCA**2)
  1023. U=U-SCA*XVEC
  1024. V=V-SCA*YVEC
  1025. W=W-SCA*ZVEC
  1026. S=SQRT(U**2+V**2+W**2)
  1027. IF (S.EQ.0.) CALL ERREUR(40)
  1028. IF (IERR.NE.0) RETURN
  1029. XPJ=XP+U-U*RAYC/S
  1030. YPJ=YP+V-V*RAYC/S
  1031. ZPJ=ZP+W-W*RAYC/S
  1032. XNJ=XJV
  1033. YNJ=YJV
  1034. ZNJ=ZJV
  1035. GOTO 350
  1036. C PROJECTION SUR CONE SOMMET XIPL DIR XIV TANG ANGLE SOMM TANGI EN
  1037. C RESTANT DANS XVEC
  1038. 230 CONTINUE
  1039. U=XIPL-XP
  1040. V=YIPL-YP
  1041. W=ZIPL-ZP
  1042. SCA=U*XIV+V*YIV+W*ZIV
  1043. U=U-SCA*XIV
  1044. V=V-SCA*YIV
  1045. W=W-SCA*ZIV
  1046. XC=XP+U
  1047. YC=YP+U
  1048. ZC=ZP+U
  1049. RL2=(SCA*TANGI)**2
  1050. SCA=U*XVEC+V*YVEC+W*ZVEC
  1051. DL2=SCA**2
  1052. IF (DL2.GT.RL2) CALL ERREUR(40)
  1053. IF(IERR.NE.0) RETURN
  1054. U=U-SCA*XVEC
  1055. V=V-SCA*YVEC
  1056. W=W-SCA*ZVEC
  1057. S=SQRT(U**2+V**2+W**2)
  1058. IF (S.EQ.0.) CALL ERREUR(40)
  1059. IF (IERR.NE.0) RETURN
  1060. U1=U/S
  1061. V1=V/S
  1062. W1=W/S
  1063. VEC=SQRT(RL2-DL2)
  1064. U=U-VEC*U1
  1065. V=V-VEC*V1
  1066. W=W-VEC*W1
  1067. XPI=XP+U
  1068. YPI=YP+V
  1069. ZPI=ZP+W
  1070. U=XIPL-XPI
  1071. V=YIPL-YPI
  1072. W=ZIPL-ZPI
  1073. U1=V*ZIV-W*YIV
  1074. V1=W*XIV-U*ZIV
  1075. W1=U*YIV-V*XIV
  1076. U2=V*W1-W*V1
  1077. V2=W*U1-U*W1
  1078. W2=U*V1-V*U1
  1079. XNI=YVEC*W2-ZVEC*V2
  1080. YNI=ZVEC*U2-XVEC*W2
  1081. ZNI=XVEC*V2-YVEC*U2
  1082. S=SQRT(XNI**2+YNI**2+ZNI**2)
  1083. IF (S.EQ.0.) CALL ERREUR(40)
  1084. IF (IERR.NE.0) RETURN
  1085. XNI=XNI/S
  1086. YNI=YNI/S
  1087. ZNI=ZNI/S
  1088. GOTO 250
  1089. 330 CONTINUE
  1090. U=XJPL-XP
  1091. V=YJPL-YP
  1092. W=ZJPL-ZP
  1093. SCA=U*XJV+V*YJV+W*ZJV
  1094. U=U-SCA*XJV
  1095. V=V-SCA*YJV
  1096. W=W-SCA*ZJV
  1097. XC=XP+U
  1098. YC=YP+U
  1099. ZC=ZP+U
  1100. RL2=(SCA*TANGJ)**2
  1101. SCA=U*XVEC+V*YVEC+W*ZVEC
  1102. DL2=SCA**2
  1103. IF (DL2.GT.RL2) CALL ERREUR(40)
  1104. IF(IERR.NE.0) RETURN
  1105. U=U-SCA*XVEC
  1106. V=V-SCA*YVEC
  1107. W=W-SCA*ZVEC
  1108. S=SQRT(U**2+V**2+W**2)
  1109. IF (S.EQ.0.) CALL ERREUR(40)
  1110. IF (IERR.NE.0) RETURN
  1111. U1=U/S
  1112. V1=V/S
  1113. W1=W/S
  1114. VEC=SQRT(RL2-DL2)
  1115. U=U-VEC*U1
  1116. V=V-VEC*V1
  1117. W=W-VEC*W1
  1118. XPJ=XP+U
  1119. YPJ=YP+V
  1120. ZPJ=ZP+W
  1121. U=XJPL-XPJ
  1122. V=YJPL-YPJ
  1123. W=ZJPL-ZPJ
  1124. U1=V*ZJV-W*YJV
  1125. V1=W*XJV-U*ZJV
  1126. W1=U*YJV-V*XJV
  1127. U2=V*W1-W*V1
  1128. V2=W*U1-U*W1
  1129. W2=U*V1-V*U1
  1130. XNJ=YVEC*W2-ZVEC*V2
  1131. YNJ=ZVEC*U2-XVEC*W2
  1132. ZNJ=XVEC*V2-YVEC*U2
  1133. S=SQRT(XNJ**2+YNJ**2+ZNJ**2)
  1134. IF (S.EQ.0.) CALL ERREUR(40)
  1135. IF (IERR.NE.0) RETURN
  1136. XNJ=XNJ/S
  1137. YNJ=YNJ/S
  1138. ZNJ=ZNJ/S
  1139. GOTO 350
  1140. C PROJECTION SUR UN TORE CENTRE XIPL AXE XIAXE CENTRE SECONDAIRE XIP1
  1141. C GRAND RAYON GRAYI PETIT RAYON PRAYI EN RESTANT DANS XVEC
  1142. 240 CONTINUE
  1143. U=XP-XIPL
  1144. V=YP-YIPL
  1145. W=ZP-ZIPL
  1146. SCA=U*XIAXE+V*YIAXE+W*ZIAXE
  1147. U=U-SCA*XIAXE
  1148. V=V-SCA*YIAXE
  1149. W=W-SCA*ZIAXE
  1150. S=SQRT(U**2+V**2+W**2)
  1151. IF (S.EQ.0.) CALL ERREUR(40)
  1152. IF (IERR.NE.0) RETURN
  1153. C CENTRE DU CERCLE INTERS
  1154. XC=XIPL+U*GRAYI/S
  1155. YC=YIPL+V*GRAYI/S
  1156. ZC=ZIPL+W*GRAYI/S
  1157. U=XP-XC
  1158. V=YP-YC
  1159. W=ZP-ZC
  1160. S=SQRT(U**2+V**2+W**2)
  1161. IF (S.EQ.0) CALL ERREUR(40)
  1162. IF (IERR.NE.0) RETURN
  1163. C PROJECTION SUR LE CERCLE
  1164. X1=XC+U*PRAYI/S
  1165. Y1=YC+V*PRAYI/S
  1166. Z1=ZC+W*PRAYI/S
  1167. C IL FAUT TOURNER SUR LE TORE JUSQU'AU PLAN XP XVEC
  1168. C CENTRE DE ROTATION
  1169. U=X1-XIPL
  1170. V=Y1-YIPL
  1171. W=Z1-ZIPL
  1172. SCA=U*XIAXE+V*YIAXE+W*ZIAXE
  1173. U=SCA*XIAXE
  1174. V=SCA*YIAXE
  1175. W=SCA*ZIAXE
  1176. XC1=XIPL+U
  1177. YC1=YIPL+V
  1178. ZC1=ZIPL+W
  1179. C RAYON DE CE CERCLE
  1180. RAY=SQRT((X1-XC1)**2+(Y1-YC1)**2+(Z1-ZC1)**2)
  1181. C DIRECTION DE L'INTERSECTION DU PLAN DU CERCLE ET DE XP XVEC
  1182. U=YIAXE*ZVEC-ZIAXE*YVEC
  1183. V=ZIAXE*XVEC-XIAXE*ZVEC
  1184. W=XIAXE*YVEC-YIAXE*XVEC
  1185. S=SQRT(U**2+V**2+W**2)
  1186. IF (S.EQ.0.) CALL ERREUR(40)
  1187. IF (IERR.NE.0) RETURN
  1188. U=U/S
  1189. V=V/S
  1190. W=W/S
  1191. U1=V*ZIAXE-W*YIAXE
  1192. V1=W*XIAXE-U*ZIAXE
  1193. W1=U*YIAXE-V*XIAXE
  1194. U2=XC1-XP
  1195. V2=YC1-YP
  1196. W2=ZC1-ZP
  1197. DNUM=U2*XVEC+V2*YVEC+W2*ZVEC
  1198. DNOM=U1*XVEC+V1*YVEC+W1*ZVEC
  1199. IF (DNOM.EQ.0.) CALL ERREUR(40)
  1200. IF (IERR.NE.0.) RETURN
  1201. PAR=-DNUM/DNOM
  1202. XC2=XC1+PAR*U1
  1203. YC2=YC1+PAR*V1
  1204. ZC2=ZC1+PAR*W1
  1205. S=SQRT((XC2-XC1)**2+(YC2-YC1)**2+(ZC2-ZC1)**2)
  1206. IF (S.GT.RAY) CALL ERREUR(40)
  1207. IF (IERR.NE.0) RETURN
  1208. D=SQRT(RAY**2-S**2)
  1209. IF ((X1-XC2)*U+(Y1-YC2)*V+(Z1-ZC2)*W.LE.0.) D=-D
  1210. XPI=XC2+D*U
  1211. YPI=YC2+D*V
  1212. ZPI=ZC2+D*W
  1213. C CENTRE DU PETIT CERCLE
  1214. U=XPI-XIPL
  1215. V=YPI-YIPL
  1216. W=ZPI-ZIPL
  1217. SCA=U*XIAXE+V*YIAXE+W*ZIAXE
  1218. U=U-SCA*XIAXE
  1219. V=V-SCA*YIAXE
  1220. W=W-SCA*ZIAXE
  1221. S=SQRT(U**2+V**2+W**2)
  1222. XC1=XIPL+U*GRAYI/S
  1223. YC1=YIPL+V*GRAYI/S
  1224. ZC1=ZIPL+W*GRAYI/S
  1225. U=XPI-XC1
  1226. V=YPI-YC1
  1227. W=ZPI-ZC1
  1228. XNI=V*ZVEC-W*YVEC
  1229. YNI=W*XVEC-U*ZVEC
  1230. ZNI=U*YVEC-V*XVEC
  1231. S=SQRT(XNI**2+YNI**2+ZNI**2)
  1232. IF (S.EQ.0.) GOTO 250
  1233. XNI=XNI/S
  1234. YNI=YNI/S
  1235. ZNI=ZNI/S
  1236. GOTO 250
  1237. 340 CONTINUE
  1238. U=XP-XJPL
  1239. V=YP-YJPL
  1240. W=ZP-ZJPL
  1241. SCA=U*XJAXE+V*YJAXE+W*ZJAXE
  1242. U=U-SCA*XJAXE
  1243. V=V-SCA*YJAXE
  1244. W=W-SCA*ZJAXE
  1245. S=SQRT(U**2+V**2+W**2)
  1246. IF (S.EQ.0.) CALL ERREUR(40)
  1247. IF (IERR.NE.0) RETURN
  1248. C CENTRE DU CERCLE INTERS
  1249. XC=XJPL+U*GRAYJ/S
  1250. YC=YJPL+V*GRAYJ/S
  1251. ZC=ZJPL+W*GRAYJ/S
  1252. U=XP-XC
  1253. V=YP-YC
  1254. W=ZP-ZC
  1255. S=SQRT(U**2+V**2+W**2)
  1256. IF (S.EQ.0) CALL ERREUR(40)
  1257. IF (IERR.NE.0) RETURN
  1258. C PROJECTION SUR LE CERCLE
  1259. X1=XC+U*PRAYJ/S
  1260. Y1=YC+V*PRAYJ/S
  1261. Z1=ZC+W*PRAYJ/S
  1262. C IL FAUT TOURNER SUR LE TORE JUSQU'AU PLAN XP XVEC
  1263. C CENTRE DE ROTATION
  1264. U=X1-XJPL
  1265. V=Y1-YJPL
  1266. W=Z1-ZJPL
  1267. SCA=U*XJAXE+V*YJAXE+W*ZJAXE
  1268. U=SCA*XJAXE
  1269. V=SCA*YJAXE
  1270. W=SCA*ZJAXE
  1271. XC1=XJPL+U
  1272. YC1=YJPL+V
  1273. ZC1=ZJPL+W
  1274. C RAYON DE CE CERCLE
  1275. RAY=SQRT((X1-XC1)**2+(Y1-YC1)**2+(Z1-ZC1)**2)
  1276. C DIRECTION DE L'INTERSECTION DU PLAN DU CERCLE ET DE XP XVEC
  1277. U=YJAXE*ZVEC-ZJAXE*YVEC
  1278. V=ZJAXE*XVEC-XJAXE*ZVEC
  1279. W=XJAXE*YVEC-YJAXE*XVEC
  1280. S=SQRT(U**2+V**2+W**2)
  1281. IF (S.EQ.0.) CALL ERREUR(40)
  1282. IF (IERR.NE.0) RETURN
  1283. U=U/S
  1284. V=V/S
  1285. W=W/S
  1286. U1=V*ZJAXE-W*YJAXE
  1287. V1=W*XJAXE-U*ZJAXE
  1288. W1=U*YJAXE-V*XJAXE
  1289. U2=XC1-XP
  1290. V2=YC1-YP
  1291. W2=ZC1-ZP
  1292. DNUM=U2*XVEC+V2*YVEC+W2*ZVEC
  1293. DNOM=U1*XVEC+V1*YVEC+W1*ZVEC
  1294. IF (DNOM.EQ.0.) CALL ERREUR(40)
  1295. IF (IERR.NE.0.) RETURN
  1296. PAR=-DNUM/DNOM
  1297. XC2=XC1+PAR*U1
  1298. YC2=YC1+PAR*V1
  1299. ZC2=ZC1+PAR*W1
  1300. S=SQRT((XC2-XC1)**2+(YC2-YC1)**2+(ZC2-ZC1)**2)
  1301. IF (S.GT.RAY) CALL ERREUR(40)
  1302. IF (IERR.NE.0) RETURN
  1303. D=SQRT(RAY**2-S**2)
  1304. IF ((X1-XC2)*U+(Y1-YC2)*V+(Z1-ZC2)*W.LE.0.) D=-D
  1305. XPJ=XC2+D*U
  1306. YPJ=YC2+D*V
  1307. ZPJ=ZC2+D*W
  1308. C CENTRE DU PETIT CERCLE
  1309. U=XPJ-XJPL
  1310. V=YPJ-YJPL
  1311. W=ZPJ-ZJPL
  1312. SCA=U*XJAXE+V*YJAXE+W*ZJAXE
  1313. U=U-SCA*XJAXE
  1314. V=V-SCA*YJAXE
  1315. W=W-SCA*ZJAXE
  1316. S=SQRT(U**2+V**2+W**2)
  1317. XC1=XJPL+U*GRAYJ/S
  1318. YC1=YJPL+V*GRAYJ/S
  1319. ZC1=ZJPL+W*GRAYJ/S
  1320. U=XPJ-XC1
  1321. V=YPJ-YC1
  1322. W=ZPJ-ZC1
  1323. XNJ=V*ZVEC-W*YVEC
  1324. YNJ=W*XVEC-U*ZVEC
  1325. ZNJ=U*YVEC-V*XVEC
  1326. S=SQRT(XNJ**2+YNJ**2+ZNJ**2)
  1327. IF (S.EQ.0.) GOTO 250
  1328. XNJ=XNJ/S
  1329. YNJ=YNJ/S
  1330. ZNJ=ZNJ/S
  1331. GOTO 350
  1332. END
  1333.  
  1334.  
  1335.  
  1336.  

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