Télécharger cadrcl.eso

Retour à la liste

Numérotation des lignes :

cadrcl
  1. C CADRCL SOURCE CB215821 23/01/25 21:15:06 11573
  2. C TRACE DE DEFORMES
  3. C CALCUL DU CADRE GLOBAL
  4. C PROJECTION SUCCESSIVE DE CHAQUE DEFORME (SELON ICLE)
  5. C 1995 option DIRE P.PEGON JRC-ISPRA
  6. C
  7. C PP option DIRE
  8. SUBROUTINE CADRCL(KABCOR,LABCO2,IOEIL,XPROJ,
  9. # ICLE,XMINT,YMINT,XMAXT,YMAXT,zmint,zmaxt,cgrav,diloc,LDIRE,
  10. > axez)
  11. IMPLICIT INTEGER(I-N)
  12. DIMENSION XO(3),XP(3),XN(3),XG(3),XV(3),UI(3),UJ(3),axez(3)
  13. dimension cgrav(*)
  14. C+PP option DIRE
  15. dimension diloc(*)
  16. C+PP
  17. COMMON /CCADRC/XN,XG,XO,UI,UJ
  18. C ATTENTION
  19. SEGMENT KABCOR(0)
  20. SEGMENT KABCO2(2,0)
  21. SEGMENT LABCO2(3,0)
  22. SEGMENT ICOR2(0)
  23. SEGMENT KXPRO2(NVEC)
  24. SEGMENT SXCORD
  25. REAL XCORD(IDIM,ITE)
  26. ENDSEGMENT
  27. SEGMENT SXCOR2
  28. REAL XCOR2(IDIM,ITE)
  29. ENDSEGMENT
  30. -INC PPARAM
  31. -INC CCOPTIO
  32. -INC SMCOORD
  33. -INC CCREEL
  34. SEGMENT XPROJ(3,ITE)
  35. SEGMENT XPRO2(3,ITE)
  36. C+PP option DIRE
  37. REAL*8 SSN
  38. LOGICAL LDIRE
  39. *dbg write(6,*) 'coucou cadrcl'
  40. C+PP
  41. C ICLE = 0 RECHERCHE DES MIN MAX SUR TOUTES LES DEFORMES
  42. C CALCUL DE LA PROJECTION
  43. XMINT=1E30
  44. YMINT=XMINT
  45. ZMINT=XMINT
  46. XMAXT=-1E30
  47. YMAXT=XMAXT
  48. ZMAXT=XMAXT
  49. NCORD=KABCOR(/1)
  50. IF (IDIM.EQ.3) GOTO 50
  51.  
  52. DO 10 ICORD=1,NCORD
  53. SXCORD=KABCOR(ICORD)
  54. * write(6,*) '1 cadrcl icord sxcord',icord,sxcord
  55. ITE=XCORD(/2)
  56. DO 19 J=1,ITE
  57. XMINT=MIN(XMINT,XCORD(1,J))
  58. XMAXT=MAX(XMAXT,XCORD(1,J))
  59. YMINT=MIN(YMINT,XCORD(2,J))
  60. YMAXT=MAX(YMAXT,XCORD(2,J))
  61. ZMINT=0
  62. ZMAXT=0
  63. 19 CONTINUE
  64. 10 CONTINUE
  65.  
  66. IF (ICLE.NE.0) GOTO 1000
  67. IF (LABCO2.EQ.0) GOTO 60
  68. DO 61 ICORD=1,NCORD
  69. IF (LABCO2(3,ICORD).EQ.0) GOTO 61
  70. KABCO2=LABCO2(1,ICORD)
  71. SXCORD=KABCOR(ICORD)
  72. * write(6,*) '2 cadrcl icord sxcord',icord,sxcord
  73. ITE=XCORD(/2)
  74. NVEC=KABCO2(/2)
  75. IF (NVEC.EQ.0) GOTO 61
  76. DO 7 IVEC=1,NVEC
  77. SXCOR2=KABCO2(1,IVEC)
  78. ICOR2=KABCO2(2,IVEC)
  79. ITE=XCOR2(/2)
  80. DO 18 J=1,ITE
  81. IF (ICOR2(J).EQ.0) GOTO 18
  82. UX=XCOR2(1,J)-XCORD(1,J)
  83. UY=XCOR2(2,J)-XCORD(2,J)
  84. U1=XCOR2(1,J)-UX/3-UY/5
  85. V1=XCOR2(2,J)-UY/3+UX/5
  86. XMINT=MIN(XMINT,U1)
  87. XMAXT=MAX(XMAXT,U1)
  88. YMINT=MIN(YMINT,V1)
  89. YMAXT=MAX(YMAXT,V1)
  90. U1=XCOR2(1,J)-UX/3+UY/5
  91. V1=XCOR2(2,J)-UY/3-UX/5
  92. XMINT=MIN(XMINT,U1)
  93. XMAXT=MAX(XMAXT,U1)
  94. YMINT=MIN(YMINT,V1)
  95. YMAXT=MAX(YMAXT,V1)
  96. XMINT=MIN(XMINT,XCOR2(1,J))
  97. XMAXT=MAX(XMAXT,XCOR2(1,J))
  98. YMINT=MIN(YMINT,XCOR2(2,J))
  99. YMAXT=MAX(YMAXT,XCOR2(2,J))
  100. 18 CONTINUE
  101. 7 CONTINUE
  102. 61 CONTINUE
  103. 60 CONTINUE
  104. * pour eviter des calculs sur des quantites non significatives,
  105. * on agrandi un peu le cadre
  106. XDIFF=XMAXT-XMINT
  107. YDIFF=YMAXT-YMINT
  108. ZDIFF=ZMAXT-ZMINT
  109. XSOMM=XMAXT+XMINT
  110. YSOMM=YMAXT+YMINT
  111. ZSOMM=ZMAXT+ZMINT
  112. XDIFF=MAX(XDIFF,YDIFF/10,ZDIFF/10)
  113. YDIFF=MAX(XDIFF/10,YDIFF,ZDIFF/10)
  114. ZDIFF=MAX(XDIFF/10,YDIFF/10,ZDIFF)
  115. XMINT=(XSOMM-XDIFF)/2
  116. XMAXT=(XSOMM+XDIFF)/2
  117. YMINT=(YSOMM-YDIFF)/2
  118. YMAXT=(YSOMM+YDIFF)/2
  119. ZMINT=(ZSOMM-ZDIFF)/2
  120. ZMAXT=(ZSOMM+ZDIFF)/2
  121. * write(6,*) '1 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt
  122. IF (ICLE.NE.0) GOTO 1000
  123. RETURN
  124. 50 CONTINUE
  125. XMINT=1E30
  126. XMAXT=-XMINT
  127. YMINT=XMINT
  128. YMAXT=-YMINT
  129. ZMINT=XMINT
  130. ZMAXT=-YMINT
  131. IREF=(IOEIL-1)*4
  132. SEGACT,MCOORD
  133. XO(1)=XCOOR(IREF+1)
  134. XO(2)=XCOOR(IREF+2)
  135. XO(3)=XCOOR(IREF+3)
  136. C POINT MOYEN
  137. DO 1 I=1,3
  138. XG(I)=0.D0
  139. 1 CONTINUE
  140. NPTOT=0
  141. DO 100 ICORD=1,NCORD
  142. SXCORD=KABCOR(ICORD)
  143. * write(6,*) '3 cadrcl icord sxcord',icord,sxcord
  144. ITE=XCORD(/2)
  145. NPTOT=NPTOT+ITE
  146. DO 2 I=1,ITE
  147. DO 3 J=1,3
  148. XG(J)=XG(J)+XCORD(J,I)
  149. 3 CONTINUE
  150. 2 CONTINUE
  151. 100 CONTINUE
  152. C+PP option DIRE
  153. IF (LDIRE)THEN
  154. DO J=1,3
  155. XG(J)=cgrav(J)
  156. XN(J)=XO(J)-XG(J)
  157. ENDDO
  158. ELSE
  159. C+PP
  160. DO 4 J=1,3
  161. XG(J)=XG(J)/NPTOT
  162. XN(J)=XO(J)-XG(J)
  163. cgrav(j)=xg(j)
  164. 4 CONTINUE
  165. C+PP option DIRE
  166. ENDIF
  167. C+PP
  168. C NORMALE AU PLAN
  169. SNI=MAX(ABS(XN(1)),ABS(XN(2)),ABS(XN(3)))+XSPETI
  170. SN=SQRT((XN(1)/SNI)**2+(XN(2)/SNI)**2+(XN(3)/SNI)**2)
  171. SN=SN*SNI
  172. IF (SN.EQ.0.) CALL ERREUR(21)
  173. IF (IERR.NE.0) RETURN
  174. C+PP option DIRE
  175. SSN=SN
  176. C+PP
  177. DO 5 J=1,3
  178. XN(J)=XN(J)/SN
  179. 5 CONTINUE
  180. C REPERE LOCAL SUR LE PLAN
  181. C+PP option DIRE
  182. IF (LDIRE)THEN
  183. DO J=1,3
  184. UJ(J)=diloc(J)
  185. ENDDO
  186. ELSE
  187. C+PP
  188. UJ(1)=0.D0
  189. UJ(2)=0.D0
  190. UJ(3)=1.D0
  191. C+PP option DIRE
  192. ENDIF
  193. C+PP
  194. 21 CONTINUE
  195. SV=UJ(1)*XN(1)+UJ(2)*XN(2)+UJ(3)*XN(3)
  196. DO 20 J=1,3
  197. UJ(J)=UJ(J)-SV*XN(J)
  198. 20 CONTINUE
  199. SV=UJ(1)**2+UJ(2)**2+UJ(3)**2
  200. IF (ABS(SV).LT.0.1) THEN
  201. UJ(1)=0.D0
  202. UJ(2)=1.D0
  203. UJ(3)=1.D0
  204. GOTO 21
  205. ENDIF
  206. SV=SQRT(SV)
  207. UJ(1)=UJ(1)/SV
  208. UJ(2)=UJ(2)/SV
  209. UJ(3)=UJ(3)/SV
  210. UI(1)=UJ(2)*XN(3)-UJ(3)*XN(2)
  211. UI(2)=UJ(3)*XN(1)-UJ(1)*XN(3)
  212. UI(3)=UJ(1)*XN(2)-UJ(2)*XN(1)
  213. C PROJECTION CONIQUE SUR LE PLAN
  214. DO 170 ICORD=1,NCORD
  215. SXCORD=KABCOR(ICORD)
  216. * write(6,*) '4 cadrcl icord sxcord',icord,sxcord
  217. ITE=XCORD(/2)
  218. DO 12 I=1,ITE
  219. DO 13 J=1,3
  220. XP(J)=XCORD(J,I)
  221. XV(J)=XP(J)-XO(J)
  222. 13 CONTINUE
  223. SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3)
  224. C+PP option DIRE
  225. IF ((LDIRE.AND.-SV.GE.SSN).OR.(.NOT.LDIRE))THEN
  226. C+PP
  227. SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3))
  228. $ *XN(3)
  229. DO 14 J=1,3
  230. XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J)
  231. 14 CONTINUE
  232. XPRO=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3)
  233. YPRO=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3)
  234. ZPRO=-SN
  235. XMINT=MIN(XMINT,XPRO)
  236. XMAXT=MAX(XMAXT,XPRO)
  237. YMINT=MIN(YMINT,YPRO)
  238. YMAXT=MAX(YMAXT,YPRO)
  239. ZMINT=MIN(ZMINT,ZPRO)
  240. ZMAXT=MAX(ZMAXT,ZPRO)
  241. C+PP option DIRE
  242. ENDIF
  243. C+PP
  244. 12 CONTINUE
  245. 170 CONTINUE
  246. * write(6,*) 'cadrcl en 239 ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt
  247. IF (LABCO2.EQ.0) GOTO 11
  248. DO 171 ICORD=1,NCORD
  249. IF (LABCO2(3,ICORD).EQ.0) GOTO 171
  250. KABCO2=LABCO2(1,ICORD)
  251. SXCORD=KABCOR(ICORD)
  252. * write(6,*) '5 cadrcl icord sxcord',icord,sxcord
  253. ITE=XCORD(/2)
  254. NVEC=KABCO2(/2)
  255. IF (NVEC.EQ.0) GOTO 171
  256. DO 6 IVEC=1,NVEC
  257. SXCOR2=KABCO2(1,IVEC)
  258. ICOR2=KABCO2(2,IVEC)
  259. DO 17 I=1,ITE
  260. IF(I.GT.ICOR2(/1)) goto 17
  261. IF (ICOR2(I).EQ.0) GOTO 17
  262. C+PP a faire meme ss DIRE car sinon que signifie XPRO,YPRO?
  263. DO J=1,3
  264. XP(J)=XCORD(J,I)
  265. XV(J)=XP(J)-XO(J)
  266. ENDDO
  267. SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3)
  268. C+PP option DIRE
  269. IF (LDIRE.AND.-SV.LT.SSN) THEN
  270. ICOR2(I)=0
  271. GOTO 17
  272. ENDIF
  273. C+PP a faire meme ss DIRE car sinon que signifie XPRO,YPRO?
  274. SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3))
  275. $ *XN(3)
  276. DO J=1,3
  277. XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J)
  278. ENDDO
  279. XPRO=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3)
  280. YPRO=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3)
  281. C+PP
  282. DO 15 J=1,3
  283. XP(J)=XCOR2(J,I)
  284. XV(J)=XP(J)-XO(J)
  285. 15 CONTINUE
  286. SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3)
  287. C+PP option DIRE
  288. IF (LDIRE.AND.-SV.LT.SSN) THEN
  289. ICOR2(I)=0
  290. GOTO 17
  291. ENDIF
  292. C+PP
  293. SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3))
  294. $ *XN(3)
  295. DO 16 J=1,3
  296. XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J)
  297. 16 CONTINUE
  298. XPROO=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3)
  299. YPROO=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3)
  300. XMINT=MIN(XMINT,XPROO)
  301. XMAXT=MAX(XMAXT,XPROO)
  302. YMINT=MIN(YMINT,YPROO)
  303. YMAXT=MAX(YMAXT,YPROO)
  304. UX=XPROO-XPRO
  305. UY=YPROO-YPRO
  306. U1=XPROO-UX/3-UY/5
  307. V1=YPROO-UY/3+UX/5
  308. XMINT=MIN(XMINT,U1)
  309. XMAXT=MAX(XMAXT,V1)
  310. U1=XPROO-UX/3+UY/5
  311. V1=YPROO-UY/3-UX/5
  312. XMINT=MIN(XMINT,U1)
  313. XMAXT=MAX(XMAXT,V1)
  314. 17 CONTINUE
  315. 6 CONTINUE
  316. 171 CONTINUE
  317. 11 CONTINUE
  318. * pour eviter des calculs sur des quantites non significatives,
  319. * on agrandi un peu le cadre
  320. * write(6,*) '2 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt
  321. XDIFF=XMAXT-XMINT
  322. YDIFF=YMAXT-YMINT
  323. ZDIFF=ZMAXT-ZMINT
  324. XSOMM=XMAXT+XMINT
  325. YSOMM=YMAXT+YMINT
  326. ZSOMM=ZMAXT+ZMINT
  327. * write(6,*) '4 cadrcl ',xdiff,ydiff,zdiff,xsomm,ysomm,zsomm
  328. XDIFF=MAX(XDIFF,YDIFF/10,ZDIFF/10)
  329. YDIFF=MAX(XDIFF/10,YDIFF,ZDIFF/10)
  330. ZDIFF=MAX(XDIFF/10,YDIFF/10,ZDIFF)
  331. * write(6,*) '5 cadrcl ',xdiff,ydiff,zdiff
  332. XMINT=(XSOMM-XDIFF)/2
  333. XMAXT=(XSOMM+XDIFF)/2
  334. YMINT=(YSOMM-YDIFF)/2
  335. YMAXT=(YSOMM+YDIFF)/2
  336. ZMINT=(ZSOMM-ZDIFF)/2
  337. ZMAXT=(ZSOMM+ZDIFF)/2
  338. * write(6,*) '3 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt
  339. * axez pour tourner correctement avec opengl
  340. axez(1)=0
  341. axez(2)=uj(3)
  342. axez(3)=sqrt(1-uj(3)**2)
  343. if (xn(3).lt.0) axez(3)=-axez(3)
  344. * write (6,*) ' axez ',axez(1),axez(2),axez(3)
  345. IF (ICLE.NE.0) GOTO 1000
  346. RETURN
  347. 1000 CONTINUE
  348. C ON REMPLIT XPROJ POUR LA DEFORME CONCERNEE
  349. SXCORD=KABCOR(ICLE)
  350. * write(6,*) '6 cadrcl icle sxcord',icord,sxcord
  351. ITE=XCORD(/2)
  352. IF (IDIM.EQ.2) GOTO 1100
  353. C+PP option DIRE
  354. SNI=MAX(ABS(XN(1)),ABS(XN(2)),ABS(XN(3)))+XSPETI
  355. SN=SQRT((XN(1)/SNI)**2+(XN(2)/SNI)**2+(XN(3)/SNI)**2)
  356. SN=SN*SNI
  357. SSN=0.95*SSN
  358. C+PP
  359. C PROJECTION CONIQUE SUR LE PLAN
  360. DO 612 I=1,ITE
  361. DO 613 J=1,3
  362. XP(J)=XCORD(J,I)
  363. XV(J)=XP(J)-XO(J)
  364. 613 CONTINUE
  365. SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3)
  366. SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3))*XN(3)
  367. XPROJ(3,I)=-SN
  368. DO 614 J=1,3
  369. XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J)
  370. 614 CONTINUE
  371. C+PP option DIRE
  372. IF ((LDIRE.AND.-SV.GE.SSN).OR.(.NOT.LDIRE))THEN
  373. C+PP
  374. XPROJ(1,I)=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3)
  375. XPROJ(2,I)=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3)
  376. C+PP option DIRE
  377. ELSE
  378. XPROJ(1,I)=(XMINT+XMAXT)/2
  379. XPROJ(2,I)=(YMINT+YMAXT)/2
  380. ENDIF
  381. C+PP
  382. 612 CONTINUE
  383. IF (LABCO2.EQ.0) GOTO 618
  384. IF (LABCO2(3,ICLE).EQ.0) GOTO 618
  385. KABCO2=LABCO2(1,ICLE)
  386. NVEC=KABCO2(/2)
  387. SEGINI KXPRO2
  388. LABCO2(2,ICLE)=KXPRO2
  389. IF (NVEC.EQ.0) GOTO 618
  390. DO 8 IVEC=1,NVEC
  391. SXCOR2=KABCO2(1,IVEC)
  392. ICOR2=KABCO2(2,IVEC)
  393. * ajout SG 2016/07/16 apparemment le ITE de XCOR2 n'est pas
  394. * forcement celui de XCORD en cas de coupe
  395. ITE=XCOR2(/2)
  396. SEGINI XPRO2
  397. KXPRO2(IVEC)=XPRO2
  398. DO 617 I=1,ITE
  399. IF (ICOR2(I).EQ.0) GOTO 617
  400. DO 615 J=1,3
  401. XP(J)=XCOR2(J,I)
  402. XV(J)=XP(J)-XO(J)
  403. 615 CONTINUE
  404. SV=XV(1)*XN(1)+XV(2)*XN(2)+XV(3)*XN(3)
  405. SN=(XP(1)-XG(1))*XN(1)+(XP(2)-XG(2))*XN(2)+(XP(3)-XG(3))
  406. $ *XN(3)
  407. XPRO2(3,I)=-SN
  408. DO 616 J=1,3
  409. XP(J)=XP(J)-(SN/SV)*XV(J)-XG(J)
  410. 616 CONTINUE
  411. * write(6,*) ' cadrcl on passe la 353 ',i
  412. XPRO2(1,I)=XP(1)*UI(1)+XP(2)*UI(2)+XP(3)*UI(3)
  413. XPRO2(2,I)=XP(1)*UJ(1)+XP(2)*UJ(2)+XP(3)*UJ(3)
  414. 617 CONTINUE
  415. ** SEGSUP SXCOR2
  416. 8 CONTINUE
  417. 618 CONTINUE
  418. ** SEGSUP SXCORD
  419. * write(6,*) '6 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt
  420. RETURN
  421. 1100 CONTINUE
  422. DO 1110 I=1,IDIM
  423. DO 1110 J=1,ITE
  424. XPROJ(I,J)=XCORD(I,J)
  425. 1110 CONTINUE
  426. IF (LABCO2.EQ.0) GOTO 1111
  427. IF (LABCO2(3,ICLE).EQ.0) GOTO 1111
  428. KABCO2=LABCO2(1,ICLE)
  429. NVEC=KABCO2(/2)
  430. SEGINI KXPRO2
  431. LABCO2(2,ICLE)=KXPRO2
  432. IF (NVEC.EQ.0) GOTO 1111
  433. DO 9 IVEC=1,NVEC
  434. SXCOR2=KABCO2(1,IVEC)
  435. ICOR2=KABCO2(2,IVEC)
  436. segact sxcor2
  437. ITE=XCOR2(/2)
  438. SEGINI XPRO2
  439. KXPRO2(IVEC)=XPRO2
  440. * write(6,*) ' cadrcl on passe aussi la 381',j
  441. DO 1113 J=1,ITE
  442. IF (ICOR2(J).EQ.0) GOTO 1113
  443. DO 1112 I=1,IDIM
  444. XPRO2(I,J)=XCOR2(I,J)
  445. 1112 CONTINUE
  446. * write(6,*) j,xcor2,(xcor2(i,j),i=1,idim)
  447. 1113 CONTINUE
  448. ** SEGSUP SXCOR2
  449. 9 CONTINUE
  450. 1111 CONTINUE
  451. ** SEGSUP SXCORD
  452. * axez pour tourner correctement avec opengl
  453. axez(1)=0
  454. axez(2)=uj(3)
  455. axez(3)=sqrt(1-uj(3)**2)
  456. if (xn(3).lt.0) axez(3)=-axez(3)
  457. * write (6,*) ' axez ',axez(1),axez(2),axez(3)
  458. * write(6,*) '7 cadrcl ',xmint,xmaxt,ymint,ymaxt,zmint,zmaxt
  459.  
  460. RETURN
  461. END
  462.  
  463.  
  464.  
  465.  
  466.  
  467.  
  468.  
  469.  
  470.  
  471.  
  472.  
  473.  
  474.  
  475.  

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