Sources Esope | Cast3M

moment
C MOMENT    SOURCE    CHAT      05/01/13    01:51:06     5004
C       =============================================
C       SUBROUTINE OF COMPUTATION OF THE MOMENTS
        SUBROUTINE MOMENT(FRQ,ETA,TE,NSPT,F,S, IG, XL0,XL1,XL2)
C       =============================================
C
C       CALLED BY POSPRESP
C       FROM A.P.+ P.P. 27/2/90+2/4/90
C
C
C       COMPUTATION OF STATIONARY MOMENTS OF POWER SPRECTRUM
C       OF THE RESPONSE OF AN OSCILATOR WITH Wn AND Eta
C       MODULUS OF THE TRANSFER FUNCTION USED:
C
C       IG=1    |H(iW)|**2=| (Wn**2+i*2*Eta*Wn*W)/
C                             ((Wn**2-W**2)+i*2*Eta**2*W*Wn) | ** 2
C                         = (4*Eta**2*Wn**2*W**2+Wn**4)/
C                           ((Wn**2-W**2)**2+4*Eta**2*W**2*Wn**2)
C     ( Transfer Function from Base Acceleration to Absolute Acceleration )
C
C       IG=2    |H(iW)|**2=W**2/((Wn**2-W**2)**2+4*Eta**2*W**2*Wn**2)
C     ( Transfer Function from Base Acceleration to Relative Velocity )
C
C       IG=3    |H(iW)|**2=1/((Wn**2-W**2)**2+4*Eta**2*W**2*Wn**2)
C     ( Transfer Function from Base Acceleration to Relative Displacement )
C
C
C       BASE ACCELERATION POWER SPECTRUM ASSUMED LINEAR BETWEEN
C       THE POINTS (F1,S1) and (F2,S2)
C
C           S(w)=C1*w+C2
C
C               ( ANALYTICAL INTEGRATION )
C
C       =============================================
      IMPLICIT INTEGER(I-N)
        IMPLICIT REAL*8 (A-H,O-Z)
        DIMENSION F(NSPT),S(NSPT)
C
        ETQD=ETA*ETA
        ETRM=SQRT(1-ETQD)
        UMMET=1.0D0-2*ETQD
        XL0=0.0D0
        XL1=0.0D0
        XL2=0.0D0
        DPI=8*ATAN(1.0D0)
        PI=DPI/2
        OMG=DPI*FRQ
        OMG2=OMG*OMG
        OMG3=OMG2*OMG
        OMG4=OMG2*OMG2
        OMG6=OMG4*OMG2
        DO 100 IESP=1,NSPT-1
        IF (((F(IESP+1)-F(IESP)).LT.1.0D-8).OR.((S(IESP).LT.1.0D-8)
     >  .AND.(S(IESP+1).LT.1.0D-8)))   GOTO 99
        S1=S(IESP)
        S2=S(IESP+1)
        X1=DPI*F(IESP)/OMG
        X2=DPI*F(IESP+1)/OMG
        XLMPQ=(LOG((X2*X2*X2*X2-2*UMMET*X2*X2+1)/
     >              (X1*X1*X1*X1-2*UMMET*X1*X1+1)))
        XLMA=(LOG((X2*X2-2*X2*ETRM+1)/
     >             (X2*X2+2*X2*ETRM+1)*
     >             (X1*X1+2*X1*ETRM+1)/
     >             (X1*X1-2*X1*ETRM+1)))
        XLMA0=(LOG((X2*X2+2*X2*ETRM+1)/
     >             (X2*X2-2*X2*ETRM+1)*
     >             (X1*X1-2*X1*ETRM+1)/
     >             (X1*X1+2*X1*ETRM+1)))
        TNGE=((X2*X2-X1*X1)*2*ETA*ETRM/
     >        (X1*X1*X2*X2-UMMET*(X1*X1+X2*X2)+1))
        TANG=(2*ETA*(X2-X1)*(1+X2*X1)/
     >        ((1-X1*X1)*(1-X2*X2)+4*ETQD*X1*X2))
C
        IF(TANG.LT.0.D0)THEN
                     ARC=PI
                     ELSE
                     ARC=0.D0
                     END IF
        IF(TNGE.LT.0.D0)THEN
                     ARCTM=PI
                     ELSE
                     ARCTM=0.D0
                     END IF
        ARC=ARC+ATAN(TANG)
        ARCTM=ARCTM+ATAN(TNGE)
C
        DDF=(F(IESP+1)-F(IESP))*DPI
        DDF2=(F(IESP+1)**2-F(IESP)**2)*DPI**2
        C1=(S2-S1)/DDF
        C2=S1-C1*F(IESP)*DPI
C
        ZZ0=(ARC/4/ETA+XLMA0/ETRM/8)/OMG3
        ZZ1=ARCTM/ETRM/ETA/OMG2/4
        ZZ2=(ARC/4/ETA+XLMA/ETRM/8)/OMG
        ZZ3=(XLMPQ+UMMET/ETA/ETRM*ARCTM)/4
C
        GOTO(1,2,3),IG
C
 1      AP0=(C1*(4*ETQD*OMG2*ZZ3+OMG4*ZZ1)+
     >       C2*(4*ETQD*OMG2*ZZ2+OMG4*ZZ0)                   )/DPI
        AP1=(C1*(4*ETQD*OMG2*DDF+
     >           OMG4*(1+8*ETQD*UMMET)*ZZ2-4*OMG6*ETQD*ZZ0)+
     >       C2*(4*ETQD*OMG2*ZZ3+OMG4*ZZ1)                   )/DPI
        AP2=(C1*(4*ETQD*OMG2*DDF2/2+
     >           OMG4*(1+8*ETQD*UMMET)*ZZ3-4*OMG6*ETQD*ZZ1)+
     >       C2*(4*ETQD*OMG2*DDF+
     >           OMG4*(1+8*ETQD*UMMET)*ZZ2-4*OMG6*ETQD*ZZ0)  )/DPI
        GOTO 10
C
C
 2      AP0=(C1*ZZ3+C2*ZZ2)/DPI
        AP1=(C1*(DDF+2*OMG2*UMMET*ZZ2-OMG4*ZZ0)+C2*ZZ3)/DPI
        AP2=(C1*(DDF2/2.+2*OMG2*UMMET*ZZ3-OMG4*ZZ1)+
     >       C2*(DDF+2*OMG2*UMMET*ZZ2-OMG4*ZZ0))/DPI
        GOTO 10
C
 3      AP0=(C1*ZZ1+C2*ZZ0)/DPI
        AP1=(C1*ZZ2+C2*ZZ1)/DPI
        AP2=(C1*ZZ3+C2*ZZ2)/DPI
        GOTO 10
C
 10     IF (AP0.GT.0.D0) THEN
                        XL0=XL0+AP0
                        END IF
        IF (AP1.GT.0.D0) THEN
                        XL1=XL1+AP1
                        END IF
        IF (AP2.GT.0.D0) THEN
                        XL2=XL2+AP2
                        END IF
99      CONTINUE
100     CONTINUE
        RETURN
        END