Fixed bugs in cov2mmpdfreg_intfc.f

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,7 +207,7 @@ 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) 
+      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      
@@ -320,16 +325,17 @@ C  Here the covariance of the problem would be initiated
  20   CONTINUE
       enddo
 
-!       hhhh=0.0d0
+      hhhh=0.0d0
-!       do 90 Iu=1,Nu
+      do Iu=1,Nu
-!       do 90 Iv=1,Nv
+      do Iv=1,Nv
 !       WRITE(10,300) Ulev(iu),Vlev(iv),UVdens(iu,iv)
-!       hhhh=hhhh+UVdens(iu,iv)
+        hhhh=hhhh+UVdens(iu,iv)
-! 90    continue
+      enddo
-!      if (nu.gt.1.and.nv.gt.1) then
+      enddo
-!      write(11,*) 'SumSum f_uv *du*dv='
+      if (nu.gt.1.and.nv.gt.1) then
-!     1,(Ulev(2)-Ulev(1))*(Vlev(2)-Vlev(1))*hhhh
+      VALUE = (Ulev(2)-Ulev(1))*(Vlev(2)-Vlev(1))*hhhh
-!      end if
+      print *,'SumSum f_uv *du*dv=', VALUE
+      end if
       
 C      sder=sqrt(XL4-XL2*XL2/XL0)
 C      cder=-XL2/sqrt(XL0)

pywafo/src/wafo/source/mreg/intfcmod.f

@@ -132,49 +132,58 @@ C
       RETURN
       END FUNCTION SPLE
 
-
+      SUBROUTINE COVG(XL0,XL2,XL4,COV1,COV2,COV3,COV,T,N)
-
-      SUBROUTINE COVG(XL0,XL2,XL4,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
+      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
-
+      integer :: I, J, II
-
+      REAL*8 :: TT, T0
 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)
-C
 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

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()
@@ -1197,29 +1204,8 @@ class SpecData1D(WafoData):
         
         dh = h[1]-h[0]
         uvdens *= dh*dh
+        uvdens = np.rot90(uvdens,-2)
 
-#        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)
         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()