Télécharger cadrcl.eso

Retour à la liste

Numérotation des lignes :

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

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