Sources Esope | Cast3M

psphe
C PSPHE     SOURCE    PV        20/03/24    21:21:01     10554          
C    CE SOUS-PROGRAMME RAMENNE UNE SPHERE SUR SES COORDONNEES PROPRES
C
      SUBROUTINE PSPHE(IOP,FER,XPROJ,NDEB,NUMNP,ICENT,tcval)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC SMCOORD
 
-INC PPARAM
-INC CCOPTIO
      real*8 tcval(*)
      SEGMENT/FER/(NFI(ITT),MAI(IPP),ITOUR)
      SEGMENT XPROJ(3,IMAX)
      SEGMENT XINT(3,max(MAI(ITOUR+1),mai(itour+2)))
*       tcval (1)   2     3     4    5     6     7      8    9
*      SAVE XVEC1,YVEC1,ZVEC1,XVEC2,YVEC2,ZVEC2,XGRAV,YGRAV,ZGRAV
*        tcval (10) 11   12   13
*      SAVE RAYNV,XINV,YINV,ZINV
      IF (IOP.EQ.2) GOTO 100
      IMCT=MAI(ITOUR+1)
      INCT=MAI(1)+1
      IMAX=(IMCT**2)/4+10
      CALL LIRENT(IMAX,0,IRETOU)
      IF (IRETOU.NE.0) IMAX=MAX(1,IMAX)
      NDEB=IMCT+1
      SEGINI XPROJ
      SEGINI XINT
      SEGACT MCOORD*mod
C   CENTRE DE LA SPHERE
      IREF=ICENT*4-3
      XCENT=XCOOR(IREF)
      YCENT=XCOOR(IREF+1)
      ZCENT=XCOOR(IREF+2)
C   CENTRE DE GRAVITE POUR DETERMINER LE CENTRE DE L'INVERSION
C    CALCUL DU RAYON DE LA SPHERE
      XGRAVS=0
      YGRAVS=0
      ZGRAVS=0
      RAYON=0.
      DO 1 I=INCT,IMCT
      IREF=NFI(I)*4-3
      XP=XCOOR(IREF)
      YP=XCOOR(IREF+1)
      ZP=XCOOR(IREF+2)
      XGRAVS=XGRAVS+XP
      YGRAVS=YGRAVS+YP
      ZGRAVS=ZGRAVS+ZP
      RAYON=RAYON+(XCENT-XP)**2+(YCENT-YP)**2+(ZCENT-ZP)**2
   1  CONTINUE
      XGRAVS=XGRAVS/(IMCT-INCT+1)
      YGRAVS=YGRAVS/(IMCT-INCT+1)
      ZGRAVS=ZGRAVS/(IMCT-INCT+1)
      RAYON=SQRT(RAYON/(IMCT-INCT+1))
 
C  CENTRE DE L'INVERSION
      XDIR=XCENT-XGRAVS
      YDIR=YCENT-YGRAVS
      ZDIR=ZCENT-ZGRAVS
      DDIR=SQRT(XDIR**2+YDIR**2+ZDIR**2)
      COF=RAYON/DDIR
      XINV=XCENT+COF*XDIR
      YINV=YCENT+COF*YDIR
      ZINV=ZCENT+COF*ZDIR
      tcval(11)=xinv
      tcval(12)=yinv
      tcval(13)=zinv
C  EN AVANT POUR L'INVERSION  ATTENTION AU CALCUL DE LA DENSITE
      RAYNV=4*RAYON**2
      tcval(10)=raynv
      DO 40 I=INCT,max(IMCT,mai(itour+2))
      II=NFI(I)
      IREF=II*4-3
      XV=XCOOR(IREF)-XINV
      YV=XCOOR(IREF+1)-YINV
      ZV=XCOOR(IREF+2)-ZINV
      DV=XV**2+YV**2+ZV**2
      IF (DV.EQ.0.) CALL ERREUR(21)
      IF (IERR.NE.0) RETURN
      XINT(1,I)=XV*RAYNV/DV
      XINT(2,I)=YV*RAYNV/DV
      XINT(3,I)=ZV*RAYNV/DV
  40  CONTINUE
C CENTRE DE GRAVITE DU PLAN DE PROJECTION
      XGRAV=0
      YGRAV=0
      ZGRAV=0
      DO 51 I=INCT,IMCT
      XGRAV=XGRAV+XINT(1,I)
      YGRAV=YGRAV+XINT(2,I)
      ZGRAV=ZGRAV+XINT(3,I)
  51  CONTINUE
      XGRAV=XGRAV/(IMCT-INCT+1)
      YGRAV=YGRAV/(IMCT-INCT+1)
      ZGRAV=ZGRAV/(IMCT-INCT+1)
      tcval(7)=xgrav
      tcval(8)=ygrav
      tcval(9)=zgrav
C  VECTEUR NORMAL
      XNORM=0
      YNORM=0
      ZNORM=0
      DO 52 IT=1,ITOUR
      IPR=MAI(IT+1)
      XV1=XINT(1,IPR)-XGRAV
      YV1=XINT(2,IPR)-YGRAV
      ZV1=XINT(3,IPR)-ZGRAV
      DO 52 I=MAI(IT-1+1)+1,MAI(IT+1)
      XV2=XINT(1,I)-XGRAV
      YV2=XINT(2,I)-YGRAV
      ZV2=XINT(3,I)-ZGRAV
      XNORM=XNORM+YV1*ZV2-ZV1*YV2
      YNORM=YNORM+ZV1*XV2-XV1*ZV2
      ZNORM=ZNORM+XV1*YV2-XV2*YV1
      XV1=XV2
      YV1=YV2
      ZV1=ZV2
  52  CONTINUE
      DNORM=SQRT (XNORM**2+YNORM**2+ZNORM**2)
      XNORM=XNORM/DNORM
      YNORM=YNORM/DNORM
      ZNORM=ZNORM/DNORM
C  FORMATION DU REPERE
      IF (ABS(XNORM).LT.0.1) THEN
       XVEC1=1-XNORM*XNORM
       YVEC1= -XNORM*YNORM
       ZVEC1= -XNORM*ZNORM
      ELSE
       XVEC1= -YNORM*XNORM
       YVEC1=1-YNORM*YNORM
       ZVEC1= -YNORM*ZNORM
      ENDIF
      DVEC1=XVEC1**2+YVEC1**2+ZVEC1**2
      DVEC1=SQRT(DVEC1)
      XVEC1=XVEC1/DVEC1
      YVEC1=YVEC1/DVEC1
      ZVEC1=ZVEC1/DVEC1
      XVEC2=YNORM*ZVEC1-ZNORM*YVEC1
      YVEC2=ZNORM*XVEC1-XNORM*ZVEC1
      ZVEC2=XNORM*YVEC1-YNORM*XVEC1
      tcval(1)=xvec1
      tcval(2)=yvec1
      tcval(3)=Zvec1
      tcval(4)=xvec2
      tcval(5)=Yvec2
      tcval(6)=Zvec2
C  EN AVANT POUR LA PROJECTION
      DO 53 I=INCT,max(IMCT,mai(itour+2))
      II=NFI(I)
      NFI(I)=I
       XRE=XINT(1,I)-XGRAV
       YRE=XINT(2,I)-YGRAV
       ZRE=XINT(3,I)-ZGRAV
       XPROJ(1,I)=XRE*XVEC1+YRE*YVEC1+ZRE*ZVEC1
      XPROJ(2,I)=XRE*XVEC2+YRE*YVEC2+ZRE*ZVEC2
       XTEST     =XRE*XNORM+YRE*YNORM+ZRE*ZNORM
      IF (ABS(XTEST).GT.DNORM*1E-2.and.i.le.imct) CALL ERREUR(21)
       IF (IERR.NE.0) RETURN
      XPROJ(3,I)=XCOOR(II*4)
  53  CONTINUE
C  VERIFIER LES DENSITES
      DO 54 IT=1,ITOUR
      II1=MAI(IT-1+1)+1
      II2=MAI(IT+1)
      DO 54 I=II1,II2
      IF (XPROJ(3,I).NE.0) GOTO 54
      IAP=I+1
      IF (IAP.GT.II2) IAP=II1
      XPROJ(3,I)=SQRT((XPROJ(1,I)-XPROJ(1,IAP))**2+(XPROJ(2,I)-XPROJ(2,
     # IAP))**2)
  54  CONTINUE
      SEGSUP XINT
      RETURN
 100  CONTINUE
C   ON RECONSTITUE LE MAILLAGE
      xvec1=tcval(1)
      yvec1=tcval(2)
      zvec1=tcval(3)
      xvec2=tcval(4)
      yvec2=tcval(5)
      zvec2=tcval(6)
      xgrav=tcval(7)
      ygrav=tcval(8)
      zgrav=tcval(9)
      raynv=tcval(10)
      xinv=tcval(11)
      yinv=tcval(12)
      zinv=tcval(13)
      SEGACT MCOORD*mod
      IF (NDEB.GT.NUMNP) GOTO 111
      NBPTA=nbpts
      NBPTS=NBPTA+NUMNP-NDEB+1
      SEGADJ MCOORD
      DO 110 I=NDEB,NUMNP
      XI=XPROJ(1,I)*XVEC1+XPROJ(2,I)*XVEC2+XGRAV
      YI=XPROJ(1,I)*YVEC1+XPROJ(2,I)*YVEC2+YGRAV
      ZI=XPROJ(1,I)*ZVEC1+XPROJ(2,I)*ZVEC2+ZGRAV
      DI=XI**2+YI**2+ZI**2
      XCOOR(NBPTA*(IDIM+1)+1)=XI*RAYNV/DI+XINV
      XCOOR(NBPTA*(IDIM+1)+2)=YI*RAYNV/DI+YINV
      XCOOR(NBPTA*(IDIM+1)+3)=ZI*RAYNV/DI+ZINV
      XCOOR((NBPTA+1)*(IDIM+1))=XPROJ(3,I)
      NBPTA=NBPTA+1
 110  CONTINUE
 111  CONTINUE
      SEGSUP XPROJ
      RETURN
      END