Télécharger psphe.eso

Retour à la liste

Numérotation des lignes :

  1. C PSPHE SOURCE CHAT 06/06/01 21:19:43 5450
  2. C CE SOUS-PROGRAMME RAMENNE UNE SPHERE SUR SES COORDONNEES PROPRES
  3. C
  4. SUBROUTINE PSPHE(IOP,FER,XPROJ,NDEB,NUMNP,ICENT,tcval)
  5. IMPLICIT INTEGER(I-N)
  6. IMPLICIT REAL*8 (A-H,O-Z)
  7. -INC SMCOORD
  8. -INC CCOPTIO
  9. real*8 tcval(*)
  10. SEGMENT/FER/(NFI(ITT),MAI(IPP),ITOUR)
  11. SEGMENT XPROJ(3,IMAX)
  12. SEGMENT XINT(3,max(MAI(ITOUR+1),mai(itour+2)))
  13. * tcval (1) 2 3 4 5 6 7 8 9
  14. * SAVE XVEC1,YVEC1,ZVEC1,XVEC2,YVEC2,ZVEC2,XGRAV,YGRAV,ZGRAV
  15. * tcval (10) 11 12 13
  16. * SAVE RAYNV,XINV,YINV,ZINV
  17. IF (IOP.EQ.2) GOTO 100
  18. IMCT=MAI(ITOUR+1)
  19. INCT=MAI(1)+1
  20. IMAX=(IMCT**2)/4+10
  21. CALL LIRENT(IMAX,0,IRETOU)
  22. IF (IRETOU.NE.0) IMAX=MAX(1,IMAX)
  23. NDEB=IMCT+1
  24. SEGINI XPROJ
  25. SEGINI XINT
  26. SEGACT MCOORD
  27. C CENTRE DE LA SPHERE
  28. IREF=ICENT*4-3
  29. XCENT=XCOOR(IREF)
  30. YCENT=XCOOR(IREF+1)
  31. ZCENT=XCOOR(IREF+2)
  32. C CENTRE DE GRAVITE POUR DETERMINER LE CENTRE DE L'INVERSION
  33. C CALCUL DU RAYON DE LA SPHERE
  34. XGRAVS=0
  35. YGRAVS=0
  36. ZGRAVS=0
  37. RAYON=0.
  38. DO 1 I=INCT,IMCT
  39. IREF=NFI(I)*4-3
  40. XP=XCOOR(IREF)
  41. YP=XCOOR(IREF+1)
  42. ZP=XCOOR(IREF+2)
  43. XGRAVS=XGRAVS+XP
  44. YGRAVS=YGRAVS+YP
  45. ZGRAVS=ZGRAVS+ZP
  46. RAYON=RAYON+(XCENT-XP)**2+(YCENT-YP)**2+(ZCENT-ZP)**2
  47. 1 CONTINUE
  48. XGRAVS=XGRAVS/(IMCT-INCT+1)
  49. YGRAVS=YGRAVS/(IMCT-INCT+1)
  50. ZGRAVS=ZGRAVS/(IMCT-INCT+1)
  51. RAYON=SQRT(RAYON/(IMCT-INCT+1))
  52.  
  53. C CENTRE DE L'INVERSION
  54. XDIR=XCENT-XGRAVS
  55. YDIR=YCENT-YGRAVS
  56. ZDIR=ZCENT-ZGRAVS
  57. DDIR=SQRT(XDIR**2+YDIR**2+ZDIR**2)
  58. COF=RAYON/DDIR
  59. XINV=XCENT+COF*XDIR
  60. YINV=YCENT+COF*YDIR
  61. ZINV=ZCENT+COF*ZDIR
  62. tcval(11)=xinv
  63. tcval(12)=yinv
  64. tcval(13)=zinv
  65. C EN AVANT POUR L'INVERSION ATTENTION AU CALCUL DE LA DENSITE
  66. RAYNV=4*RAYON**2
  67. tcval(10)=raynv
  68. DO 40 I=INCT,max(IMCT,mai(itour+2))
  69. II=NFI(I)
  70. IREF=II*4-3
  71. XV=XCOOR(IREF)-XINV
  72. YV=XCOOR(IREF+1)-YINV
  73. ZV=XCOOR(IREF+2)-ZINV
  74. DV=XV**2+YV**2+ZV**2
  75. IF (DV.EQ.0.) CALL ERREUR(21)
  76. IF (IERR.NE.0) RETURN
  77. XINT(1,I)=XV*RAYNV/DV
  78. XINT(2,I)=YV*RAYNV/DV
  79. XINT(3,I)=ZV*RAYNV/DV
  80. 40 CONTINUE
  81. C CENTRE DE GRAVITE DU PLAN DE PROJECTION
  82. XGRAV=0
  83. YGRAV=0
  84. ZGRAV=0
  85. DO 51 I=INCT,IMCT
  86. XGRAV=XGRAV+XINT(1,I)
  87. YGRAV=YGRAV+XINT(2,I)
  88. ZGRAV=ZGRAV+XINT(3,I)
  89. 51 CONTINUE
  90. XGRAV=XGRAV/(IMCT-INCT+1)
  91. YGRAV=YGRAV/(IMCT-INCT+1)
  92. ZGRAV=ZGRAV/(IMCT-INCT+1)
  93. tcval(7)=xgrav
  94. tcval(8)=ygrav
  95. tcval(9)=zgrav
  96. C VECTEUR NORMAL
  97. XNORM=0
  98. YNORM=0
  99. ZNORM=0
  100. DO 52 IT=1,ITOUR
  101. IPR=MAI(IT+1)
  102. XV1=XINT(1,IPR)-XGRAV
  103. YV1=XINT(2,IPR)-YGRAV
  104. ZV1=XINT(3,IPR)-ZGRAV
  105. DO 52 I=MAI(IT-1+1)+1,MAI(IT+1)
  106. XV2=XINT(1,I)-XGRAV
  107. YV2=XINT(2,I)-YGRAV
  108. ZV2=XINT(3,I)-ZGRAV
  109. XNORM=XNORM+YV1*ZV2-ZV1*YV2
  110. YNORM=YNORM+ZV1*XV2-XV1*ZV2
  111. ZNORM=ZNORM+XV1*YV2-XV2*YV1
  112. XV1=XV2
  113. YV1=YV2
  114. ZV1=ZV2
  115. 52 CONTINUE
  116. DNORM=SQRT (XNORM**2+YNORM**2+ZNORM**2)
  117. XNORM=XNORM/DNORM
  118. YNORM=YNORM/DNORM
  119. ZNORM=ZNORM/DNORM
  120. C FORMATION DU REPERE
  121. IF (ABS(XNORM).LT.0.1) THEN
  122. XVEC1=1-XNORM*XNORM
  123. YVEC1= -XNORM*YNORM
  124. ZVEC1= -XNORM*ZNORM
  125. ELSE
  126. XVEC1= -YNORM*XNORM
  127. YVEC1=1-YNORM*YNORM
  128. ZVEC1= -YNORM*ZNORM
  129. ENDIF
  130. DVEC1=XVEC1**2+YVEC1**2+ZVEC1**2
  131. DVEC1=SQRT(DVEC1)
  132. XVEC1=XVEC1/DVEC1
  133. YVEC1=YVEC1/DVEC1
  134. ZVEC1=ZVEC1/DVEC1
  135. XVEC2=YNORM*ZVEC1-ZNORM*YVEC1
  136. YVEC2=ZNORM*XVEC1-XNORM*ZVEC1
  137. ZVEC2=XNORM*YVEC1-YNORM*XVEC1
  138. tcval(1)=xvec1
  139. tcval(2)=yvec1
  140. tcval(3)=Zvec1
  141. tcval(4)=xvec2
  142. tcval(5)=Yvec2
  143. tcval(6)=Zvec2
  144. C EN AVANT POUR LA PROJECTION
  145. DO 53 I=INCT,max(IMCT,mai(itour+2))
  146. II=NFI(I)
  147. NFI(I)=I
  148. XRE=XINT(1,I)-XGRAV
  149. YRE=XINT(2,I)-YGRAV
  150. ZRE=XINT(3,I)-ZGRAV
  151. XPROJ(1,I)=XRE*XVEC1+YRE*YVEC1+ZRE*ZVEC1
  152. XPROJ(2,I)=XRE*XVEC2+YRE*YVEC2+ZRE*ZVEC2
  153. XTEST =XRE*XNORM+YRE*YNORM+ZRE*ZNORM
  154. IF (ABS(XTEST).GT.DNORM*1E-2.and.i.le.imct) CALL ERREUR(21)
  155. IF (IERR.NE.0) RETURN
  156. XPROJ(3,I)=XCOOR(II*4)
  157. 53 CONTINUE
  158. C VERIFIER LES DENSITES
  159. DO 54 IT=1,ITOUR
  160. II1=MAI(IT-1+1)+1
  161. II2=MAI(IT+1)
  162. DO 54 I=II1,II2
  163. IF (XPROJ(3,I).NE.0) GOTO 54
  164. IAP=I+1
  165. IF (IAP.GT.II2) IAP=II1
  166. XPROJ(3,I)=SQRT((XPROJ(1,I)-XPROJ(1,IAP))**2+(XPROJ(2,I)-XPROJ(2,
  167. # IAP))**2)
  168. 54 CONTINUE
  169. SEGSUP XINT
  170. RETURN
  171. 100 CONTINUE
  172. C ON RECONSTITUE LE MAILLAGE
  173. xvec1=tcval(1)
  174. yvec1=tcval(2)
  175. zvec1=tcval(3)
  176. xvec2=tcval(4)
  177. yvec2=tcval(5)
  178. zvec2=tcval(6)
  179. xgrav=tcval(7)
  180. ygrav=tcval(8)
  181. zgrav=tcval(9)
  182. raynv=tcval(10)
  183. xinv=tcval(11)
  184. yinv=tcval(12)
  185. zinv=tcval(13)
  186. SEGACT MCOORD
  187. IF (NDEB.GT.NUMNP) GOTO 111
  188. NBPTA=XCOOR(/1)/(IDIM+1)
  189. NBPTS=NBPTA+NUMNP-NDEB+1
  190. SEGADJ MCOORD
  191. DO 110 I=NDEB,NUMNP
  192. XI=XPROJ(1,I)*XVEC1+XPROJ(2,I)*XVEC2+XGRAV
  193. YI=XPROJ(1,I)*YVEC1+XPROJ(2,I)*YVEC2+YGRAV
  194. ZI=XPROJ(1,I)*ZVEC1+XPROJ(2,I)*ZVEC2+ZGRAV
  195. DI=XI**2+YI**2+ZI**2
  196. XCOOR(NBPTA*(IDIM+1)+1)=XI*RAYNV/DI+XINV
  197. XCOOR(NBPTA*(IDIM+1)+2)=YI*RAYNV/DI+YINV
  198. XCOOR(NBPTA*(IDIM+1)+3)=ZI*RAYNV/DI+ZINV
  199. XCOOR((NBPTA+1)*(IDIM+1))=XPROJ(3,I)
  200. NBPTA=NBPTA+1
  201. 110 CONTINUE
  202. 111 CONTINUE
  203. SEGSUP XPROJ
  204. RETURN
  205. END
  206.  
  207.  
  208.  
  209.  

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