diff --git a/pywafo/src/wafo/cov2mod.pyd b/pywafo/src/wafo/cov2mod.pyd index edd44c3..d65754f 100644 Binary files a/pywafo/src/wafo/cov2mod.pyd and b/pywafo/src/wafo/cov2mod.pyd differ diff --git a/pywafo/src/wafo/source/mreg/cov2mmpdfreg_intfc.f b/pywafo/src/wafo/source/mreg/cov2mmpdfreg_intfc.f index 94b7e8b..3a8f89e 100644 --- a/pywafo/src/wafo/source/mreg/cov2mmpdfreg_intfc.f +++ b/pywafo/src/wafo/source/mreg/cov2mmpdfreg_intfc.f @@ -65,7 +65,7 @@ Cf2py real*8, intent(out), depend(Nu,Nv) :: UVdens Cf2py depend(Ng) Xg Cf2py depend(Nt,5) COV real*8 Q0,SQ0,Q1,SQ1, U,V,VV, XL0, XL2, XL4 - REAL*8 VDERI, CDER,SDER, DER, CONST1, F, HHHH,FM, VALUE + REAL*8 VDERI, CDER,SDER, DER, CONST1, F, HHHH, FM, VALUE C INTEGER, PARAMETER :: MMAX = 5, NMAX = 101, RDIM = 10201 REAL*8, DIMENSION(NMAX) :: HHT,VT,UT,Vdd,Udd REAL*8, DIMENSION(RDIM) :: R,R1,R2,R3 @@ -118,7 +118,7 @@ c c OBS. we are using the variables R,R1,R2 R3 as a temporary storage C for transformation g of the process. -c + N = Nt CALL INITLEVELS(T,HHT,Nt,NU,Nv) C CALL INITLEVELS(Ulev,NU,Vlev,NV,T,HHT,Nt,R1,R2,NG) IF( Tg(1) .gt. Tg(ng)) then @@ -127,8 +127,13 @@ C CALL INITLEVELS(Ulev,NU,Vlev,NV,T,HHT,Nt,R1,R2,NG) end if if(abs(Tg(ng)-Tg(1))*abs(Xg(ng)-Xg(1)).lt.0.01d0) then print *,'The transformation g is singular, stop' - stop + return end if + +! do IV=1,Nt +! print *, 'Cov', COV(IV,:) +! end do + DO IV=1,Nv V=Vlev(IV) CALL TRANSF(NG,V,Xg,Tg,VALUE,DER) @@ -145,19 +150,19 @@ C CALL INITLEVELS(Ulev,NU,Vlev,NV,T,HHT,Nt,R1,R2,NG) enddo enddo - CALL COVG(XL0,XL2,XL4,COV,T,Nt) - + CALL COVG(XL0,XL2,XL4,R1,R2,R3,COV,T,Nt) + Q0=XL4 IF (Q0.le.1.0D0+EPS) then Print *,'Covariance structure is singular, stop.' - stop + return end if SQ0 = SQRT(Q0) Q1 = XL0-XL2*XL2/XL4 IF (Q1.le.EPS) then Print *,'Covariance structure is singular, stop.' - stop + return end if SQ1 = SQRT(Q1) DO I=1,Nt @@ -202,8 +207,8 @@ c10 CONTINUE C C R1 contains Cov(X(T(I)),X'(T(J))|X'(0),X''(0),X(0)) C - R1(J+(I-1)*N)=R1(J+(I-1)*N) - COV(I,2)*(COV(J,3)/XL2) - 1 - (B0(I)*DB0(J)/Q0) - (B1(I)*DB1(J)/Q1) + R1(J+(I-1)*N) = R1(J+(I-1)*N) - COV(I,2)*(COV(J,3)/XL2) + 1 - (B0(I)*DB0(J)/Q0) - (B1(I)*DB1(J)/Q1) C C R2 contains Cov(X'(T(I)),X'(T(J))|X'(0),X''(0),X(0)) @@ -320,16 +325,17 @@ C Here the covariance of the problem would be initiated 20 CONTINUE enddo -! hhhh=0.0d0 -! do 90 Iu=1,Nu -! do 90 Iv=1,Nv + hhhh=0.0d0 + do Iu=1,Nu + do Iv=1,Nv ! WRITE(10,300) Ulev(iu),Vlev(iv),UVdens(iu,iv) -! hhhh=hhhh+UVdens(iu,iv) -! 90 continue -! if (nu.gt.1.and.nv.gt.1) then -! write(11,*) 'SumSum f_uv *du*dv=' -! 1,(Ulev(2)-Ulev(1))*(Vlev(2)-Vlev(1))*hhhh -! end if + hhhh=hhhh+UVdens(iu,iv) + enddo + enddo + if (nu.gt.1.and.nv.gt.1) then + VALUE = (Ulev(2)-Ulev(1))*(Vlev(2)-Vlev(1))*hhhh + print *,'SumSum f_uv *du*dv=', VALUE + end if C sder=sqrt(XL4-XL2*XL2/XL0) C cder=-XL2/sqrt(XL0) diff --git a/pywafo/src/wafo/source/mreg/intfcmod.f b/pywafo/src/wafo/source/mreg/intfcmod.f index bf90003..9cda573 100644 --- a/pywafo/src/wafo/source/mreg/intfcmod.f +++ b/pywafo/src/wafo/source/mreg/intfcmod.f @@ -132,49 +132,58 @@ C RETURN END FUNCTION SPLE - - - SUBROUTINE COVG(XL0,XL2,XL4,COV,T,N) + SUBROUTINE COVG(XL0,XL2,XL4,COV1,COV2,COV3,COV,T,N) C -C COVG evaluates: +C Covariance function and its four derivatives for a vector T of length N +C is assumed in a vector COV; COV(1,...,N,1)=r(T), COV(1,...,N, 2)=r'(T), etc. +C The vector COV should be of the shape N x 5. C +C COVG Returns: C XL0,XL2,XL4 - spectral moments. C -C Covariance function and its four derivatives for a vector T of length N. -C It is saved in a vector COV; COV(1,...,N)=r(T), COV(N+1,...,2N)=r'(T), etc. -C The vector COV should be of the length 5*N. -C C Covariance matrices COV1=r'(T-T), COV2=r''(T-T) and COV3=r'''(T-T) -C Dimension of COV1, COV2 should be N*N. +C Dimension of COV1, COV2 should be atleast N*N. C -! USE SIZEMOD -! IMPLICIT NONE -C INTEGER, PARAMETER:: NMAX = 101, RDIM = 10201 + USE SIZEMOD +! IMPLICIT NONE +C INTEGER, PARAMETER:: NMAX = 101, RDIM = 10201 REAL*8, PARAMETER:: ZERO = 0.0d0 REAL*8, intent(inout) :: XL0,XL2,XL4 REAL*8, DIMENSION(N,5), intent(in) :: COV REAL*8, DIMENSION(N), intent(in) :: T + REAL*8, DIMENSION(RDIM), intent(inout) :: COV1,COV2,COV3 INTEGER, intent(in) :: N - - -C -C COV(Y(T),Y(0)) = COV(:,1) -C - XL0 = COV(1,1) -! XL0 = SPLE(NT,ZERO,COV(:,1),T) -C -C DERIVATIVE COV(Y(T),Y(0)) = COV(:,2) + integer :: I, J, II + REAL*8 :: TT, T0 C +C COV(Y(T),Y(0)) = COV(:,1) +C DERIVATIVE COV(Y(T),Y(0)) = COV(:,2) C 2-DERIVATIVE COV(Y(T),Y(0)) = COV(:,3) - XL2 = -COV(1,3) -! XL2 = -SPLE(NT,ZERO,COV(:,3),T) C 3-DERIVATIVE COV(Y(T),Y(0)) = COV(:,4) - C 4-DERIVATIVE COV(Y(T),Y(0)) = COV(:,5) + XL0 = COV(1,1) + XL2 = -COV(1,3) XL4 = COV(1,5) -! XL4 = SPLE(NT,ZERO,COV(:,5),T) - +! XL0 = SPLE(NT, ZERO, COV(:,1), T) +! XL2 = -SPLE(NT, ZERO, COV(:,3), T) +! XL4 = SPLE(NT, ZERO, COV(:,5), T) + + II=0 + DO I=1,N + DO J=1,N + II = II+1 + T0 = T(J)-T(I) + TT = ABS(T0) + COV1(II) = SPLE(N, TT, COV(:,2), T) + COV2(II) = SPLE(N, TT, COV(:,3), T) + COV3(II) = SPLE(N, TT, COV(:,4), T) + IF (T0.LT.0.0d0) then + COV1(II)=-COV1(II) + COV3(II)=-COV3(II) + endif + enddo + enddo RETURN END SUBROUTINE COVG END module intfcmod diff --git a/pywafo/src/wafo/spectrum/core.py b/pywafo/src/wafo/spectrum/core.py index 38ba7af..79a111d 100644 --- a/pywafo/src/wafo/spectrum/core.py +++ b/pywafo/src/wafo/spectrum/core.py @@ -1148,6 +1148,13 @@ class SpecData1D(WafoData): NIT=5, leads tp call RIND5, IAC is maximally 4-dimensional integral, while IAC=1 leads to maximum 5 dimensional integral. + >>> import numpy as np + >>> import wafo.spectrum.models as sm + >>> Sj = sm.Jonswap(Hm0=3) + >>> w = np.linspace(0,4,256) + >>> S1 = Sj.tospecdata(w) #Make spectrum object from numerical values + >>> S = sm.SpecData1D(Sj(w),w) # Alternatively do it manually + S.to_mm_pdf() ''' S = self.copy() @@ -1177,7 +1184,7 @@ class SpecData1D(WafoData): unused_t0, tn, Nt = paramt - t = linspace(0, tn/A, Nt) # normalized times + t = linspace(0, tn/A, Nt) # normalized times #Transform amplitudes to Gaussian levels: h = linspace(*paramu); @@ -1197,29 +1204,8 @@ class SpecData1D(WafoData): dh = h[1]-h[0] uvdens *= dh*dh - -# if (defnr==0) -# f.f =fliplr(mctp2rfc(fliplr(ftmp)));%* sqrt(-R(1,6)/R(1,4))/2/pi; -# f.title ='Joint density of maximum and rainflow minimum'; -# f.labx{1}='max [m]'; -# f.labx{2}='rainflow min [m]'; -# elseif (defnr==-1) -# %CC= normalizing constant= 1/ expected number of u-up-crossings of X -# %CC = 2*pi*sqrt(L0/L2)*exp(0.5D0*u*u/L0); -# % CC = normalizing constant = 1/ expected number of zero-up-crossings of X' -# %CC = 2*pi*sqrt(L2/L4); -# fact = sqrt(L0/L4); -# -# f.f = fliplr(mctp2tc(fliplr(ftmp*fact),utc,paramu)); -# index1 = find(f.x{1}>0); -# index2 = find(f.x{2}<0); -# f.f = flipud(f.f(index2,index1)); -# f.x{1} = f.x{1}(index1); -# f.x{2} = abs(flipud(f.x{2}(index2))); -# f.title ='Joint density of crest and trough'; -# f.labx{1}='Crest [m]'; -# f.labx{2}='Trough [m]'; -# else %(defnr==1) + uvdens = np.rot90(uvdens,-2) + mmpdf = WafoData(uvdens,args=(h,h), title='Joint density of maximum and minimum', xlab='max [m]',ylab='min [m]') return mmpdf @@ -3287,9 +3273,19 @@ def main(): pylab.close('all') print('done') +def test_mm_pdf(): + import numpy as np + import wafo.spectrum.models as sm + Sj = sm.Jonswap(Hm0=7, Tp=11) + w = np.linspace(0,4,256) + S1 = Sj.tospecdata(w) #Make spectrum object from numerical values + S = sm.SpecData1D(Sj(w),w) # Alternatively do it manually + mm = S.to_mm_pdf() + if __name__ == '__main__': - if True: #False : # + if False : #True: # import doctest doctest.testmod() else: - main() + test_mm_pdf() + #main()