Télécharger cadrcl.eso

Retour à la liste

Numérotation des lignes :

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

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