diff --git a/wafo/source/mreg/cov2mmpdfreg_intfc.f b/wafo/source/mreg/cov2mmpdfreg_intfc.f index 3a8f89e..445d617 100644 --- a/wafo/source/mreg/cov2mmpdfreg_intfc.f +++ b/wafo/source/mreg/cov2mmpdfreg_intfc.f @@ -1,8 +1,8 @@ -C Version 1994-X-18 +C Version 1994-X-18 -C This is a new version of WAMP program computing crest-trough wavelength +C This is a new version of WAMP program computing crest-trough wavelength C and amplitude density. -C +C C revised pab 2007 C -moved all common blocks into modules C -renamed from minmax to sp2mmpdfreg + fixed some bugs @@ -64,8 +64,9 @@ Cf2py integer, optional :: NIT = 2 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 Q0,SQ0,Q1,SQ1, U,V, XL0, XL2, XL4 + REAL*8 VDERI, DER, F, HHHH, VALUE +C REAL*8 VV, CDER,SDER, CONST1, FM 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 @@ -79,14 +80,14 @@ C DIMENSION COV(5*NMAX),R(RDIM),R1(RDIM),R2(RDIM),R3(RDIM) C C The program computes the joint density of maximum the following minimum -C and the distance between Max and min for a zero-mean stationary -C Gaussian process with covariance function defined explicitely with 4 -C derivatives. The process should be normalized so that the first and +C and the distance between Max and min for a zero-mean stationary +C Gaussian process with covariance function defined explicitely with 4 +C derivatives. The process should be normalized so that the first and C the second spectral moments are equal to 1. The values of Max are taken C as the nodes at Hermite-Quadrature and then integrated out so that C the output is a joint density of wavelength T and amplitude H=Max-min. C The Max values are defined by subroutine Gauss_M with the accuracy -C input epsu. The principle is that the integral of the marginal density +C input epsu. The principle is that the integral of the marginal density C of f_Max is computed with sufficient accuracy. C REAL*8, DIMENSION(NMAX) :: B0,DB0,DDB0,B1,DB1,DDB1,DB2,DDB2 @@ -96,13 +97,13 @@ C DIMENSION B1(NMAX),DB1(NMAX),DDB1(NMAX) C DIMENSION DB2(NMAX),DDB2(NMAX) C DIMENSION Q(NMAX),SQ(NMAX),VDER(NMAX),DBI(NMAX),BI(NMAX) INTEGER :: J,I,I1,I2,I3,IU, IV,N, NNIT, INF - INTEGER :: fffff +C INTEGER :: fffff C REAL*8 EPS0 C INTEGER III01,III11,III21,III31,III41,III51 C *,III61,III71,III81,III91,III101 , III0 C COMMON/CHECK1/III01,III11,III21,III31,III41,III51 -C *,III61,III71,III81,III91,III101 -C COMMON/CHECKQ/III0 +C *,III61,III71,III81,III91,III101 +C COMMON/CHECKQ/III0 C COMMON /EPS/ EPS,EPSS,CEPSS C @@ -115,7 +116,7 @@ C ! OPEN(UNIT=10,FILE='Maxmin.out') ! OPEN(UNIT=11,FILE='Maxmin.log') c -c OBS. we are using the variables R,R1,R2 R3 as a temporary storage +c OBS. we are using the variables R,R1,R2 R3 as a temporary storage C for transformation g of the process. N = Nt @@ -136,22 +137,22 @@ C CALL INITLEVELS(Ulev,NU,Vlev,NV,T,HHT,Nt,R1,R2,NG) DO IV=1,Nv V=Vlev(IV) - CALL TRANSF(NG,V,Xg,Tg,VALUE,DER) + CALL TRANSF(NG,V,Xg,Tg,VALUE,DER) VT(IV)=VALUE Vdd(IV)=DER enddo DO IU=1,Nu U = Ulev(IU) - CALL TRANSF(NG,U,Xg,Tg,VALUE,DER) + CALL TRANSF(NG,U,Xg,Tg,VALUE,DER) UT(IU) = VALUE Udd(IU) = DER do IV=1,Nv UVdens(IU,IV)=0.0d0 enddo enddo - + CALL COVG(XL0,XL2,XL4,R1,R2,R3,COV,T,Nt) - + Q0=XL4 IF (Q0.le.1.0D0+EPS) then @@ -169,7 +170,7 @@ C CALL INITLEVELS(Ulev,NU,Vlev,NV,T,HHT,Nt,R1,R2,NG) B0(I) =-COV(I,3) DB0(I) =-COV(I,4) DDB0(I)=-COV(I,5) - + B1(I) =COV(I,1)+COV(I,3)*(XL2/XL4) DB1(I) =COV(I,2)+COV(I,4)*(XL2/XL4) DDB1(I)=COV(I,3)+XL2*(COV(I,5)/XL4) @@ -180,8 +181,8 @@ C Q(I)=XL0 - COV(I,2)*(COV(I,2)/XL2) - B0(I)*(B0(I)/Q0) 1 -B1(I)*(B1(I)/Q1) VDER(I)=XL4 - (COV(I,4)*COV(I,4))/XL2 - (DDB0(I)*DDB0(I))/Q0 - 1 - (DDB1(I)*DDB1(I))/Q1 - + 1 - (DDB1(I)*DDB1(I))/Q1 + C C DDB2(I) contains Cov(X''(T(i)),X(T(i))|X'(0),X''(0),X(0)) @@ -193,84 +194,84 @@ C DDB2(i)=0.0d0 else SQ(I)=SQRT(Q(I)) -C +C C VDER(I) contains Var(X''(T(i))|X'(0),X''(0),X(0),X(T(i)) C - + VDER(I)=VDER(I) - (DDB2(I)*DDB2(I))/Q(I) end if - + c10 CONTINUE enddo DO I=1,Nt DO J=1,Nt -C -C R1 contains Cov(X(T(I)),X'(T(J))|X'(0),X''(0),X(0)) +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) - -C + +C C R2 contains Cov(X'(T(I)),X'(T(J))|X'(0),X''(0),X(0)) C - R2(J+(I-1)*N) = -R2(J+(I-1)*N) - COV(I,3)*(COV(J,3)/XL2) - 1 - DB0(I)*DB0(J)/Q0 - DB1(I)*(DB1(J)/Q1) -C + R2(J+(I-1)*N) = -R2(J+(I-1)*N) - COV(I,3)*(COV(J,3)/XL2) + 1 - DB0(I)*DB0(J)/Q0 - DB1(I)*(DB1(J)/Q1) +C C R3 contains Cov(X''(T(I)),X'(T(J))|X'(0),X''(0),X(0)) C - R3(J+(I-1)*N) = R3(J+(I-1)*N) - COV(I,4)*(COV(J,3)/XL2) - 1 - DB0(J)*(DDB0(I)/Q0) - DDB1(I)*(DB1(J)/Q1) + R3(J+(I-1)*N) = R3(J+(I-1)*N) - COV(I,4)*(COV(J,3)/XL2) + 1 - DB0(J)*(DDB0(I)/Q0) - DDB1(I)*(DB1(J)/Q1) c15 CONTINUE enddo enddo C The initiations are finished and we are beginning with 3 loops C on T=T(I), U=Ulevels(IU), V=Ulevels(IV), U>V. - + DO I=1,Nt NNIT=NIT IF (Q(I).LE.EPS) GO TO 20 DO I1=1,I - DB2(I1)=R1(I1+(I-1)*N) + DB2(I1)=R1(I1+(I-1)*N) -C Cov(X'(T(I1)),X(T(i))|X'(0),X''(0),X(0)) +C Cov(X'(T(I1)),X(T(i))|X'(0),X''(0),X(0)) C DDB2(I) contains Cov(X''(T(i)),X(T(i))|X'(0),X''(0),X(0)) enddo - + DO I3=1,I - DBI(I3) = R3(I3+(I-1)*N) - (DDB2(I)*DB2(I3)/Q(I)) + DBI(I3) = R3(I3+(I-1)*N) - (DDB2(I)*DB2(I3)/Q(I)) BI(I3) = R2(I3+(I-1)*N) - (DB2(I)*DB2(I3)/Q(I)) enddo DO I3=1,I-1 AI(I3)=0.0d0 - AI(I3+I-1)=DB0(I3)/SQ0 + AI(I3+I-1)=DB0(I3)/SQ0 AI(I3+2*(I-1))=DB1(I3)/SQ1 AI(I3+3*(I-1))=DB2(I3)/SQ(I) enddo VDERI=VDER(I) DAI(1)=0.0d0 DAI(2)=DDB0(I)/SQ0 - DAI(3)=DDB1(I)/SQ1 - DAI(4)=DDB2(I)/SQ(I) - AA(1,1)=DB0(I)/SQ0 - AA(1,2)=DB1(I)/SQ1 - AA(1,3)=DB2(I)/SQ(I) + DAI(3)=DDB1(I)/SQ1 + DAI(4)=DDB2(I)/SQ(I) + AA(1,1)=DB0(I)/SQ0 + AA(1,2)=DB1(I)/SQ1 + AA(1,3)=DB2(I)/SQ(I) AA(2,1)=XL2/SQ0 AA(2,2)=SQ1 AA(2,3)=0.0d0 - AA(3,1)=B0(I)/SQ0 - AA(3,2)=B1(I)/SQ1 + AA(3,1)=B0(I)/SQ0 + AA(3,2)=B1(I)/SQ1 AA(3,3)=SQ(I) IF (BI(I).LE.EPS) NNIT=0 IF (NNIT.GT.1) THEN IF(I.LT.1) GO TO 41 DO I1=1,I-1 DO I2=1,I-1 -C R contains Cov(X'(T(I1)),X'(T(I2))|X'(0),X''(0),X(0),X(I)) - R(I2+(I1-1)*(I-1))=R2(I2+(I1-1)*N)-(DB2(I1)*DB2(I2)/Q(I)) +C R contains Cov(X'(T(I1)),X'(T(I2))|X'(0),X''(0),X(0),X(I)) + R(I2+(I1-1)*(I-1))=R2(I2+(I1-1)*N)-(DB2(I1)*DB2(I2)/Q(I)) enddo enddo @@ -283,27 +284,27 @@ C Here the covariance of the problem would be initiated Print *,' Laps to go:',Nt-I+1 DO IV=1,Nv V=VT(IV) -! IF (ABS(V).GT.5.0D0) GO TO 80 +! IF (ABS(V).GT.5.0D0) GO TO 80 IF (Vdd(IV).LT.EPS0) GO TO 80 DO IU=1,Nu U=UT(IU) IF (U.LE.V) go to 60 -! IF (ABS(U).GT.5.0D0) GO TO 60 +! IF (ABS(U).GT.5.0D0) GO TO 60 IF (Udd(IU).LT.EPS0) GO TO 60 BB(1)=0.0d0 BB(2)=U BB(3)=V -! if (IV.EQ.2.AND.IU.EQ.1) THEN +! if (IV.EQ.2.AND.IU.EQ.1) THEN ! fffff = 10 ! endif - + CALL MREG(F,R,BI,DBI,AA,BB,AI,DAI,VDERI,3,I-1,NNIT,INF) INF=1 UVdens(IU,IV) = UVdens(IU,IV) + Udd(IU)*Vdd(IV)*HHT(I)*F -! if (F.GT.0.01.AND.U.GT.2.AND.V.LT.-2) THEN +! if (F.GT.0.01.AND.U.GT.2.AND.V.LT.-2) THEN ! if (N-I+1 .eq. 38.and.IV.EQ.26.AND.IU.EQ.16) THEN -! if (IV.EQ.32.AND.IU.EQ.8.and.I.eq.11) THEN +! if (IV.EQ.32.AND.IU.EQ.8.and.I.eq.11) THEN ! PRINT * ,' R:', R(1:I) ! PRINT * ,' BI:', BI(1:I) ! PRINT * ,' DBI:', DBI(1:I) @@ -336,21 +337,21 @@ C Here the covariance of the problem would be initiated 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) C const1=1/sqrt(XL0*XL4) -C DO 95 IU=1,NU +C DO 95 IU=1,NU C U=UT(IU) C FM=Udd(IU)*const1*exp(-0.5*U*U/XL0)*PMEAN(-cder*U,sder) C WRITE(9,300) Ulev(IU),FM -C 95 continue -C DO 105 IV=1,NV +C 95 continue +C DO 105 IV=1,NV C V=VT(IV) C VV=cder*V C Fm=Vdd(IV)*const1*exp(-0.5*V*V/XL0)*PMEAN(VV,sder) C WRITE(8,300) Vlev(IV),Fm -C 105 continue +C 105 continue if (III0.eq.0) III0=1 PRINT *, 'Rate of calls RINDT0:',float(iii01)/float(III0) diff --git a/wafo/source/mvnprd/mvnprd_interface.f b/wafo/source/mvnprd/mvnprd_interface.f index 93e4b3a..ab89396 100644 --- a/wafo/source/mvnprd/mvnprd_interface.f +++ b/wafo/source/mvnprd/mvnprd_interface.f @@ -1,16 +1,17 @@ subroutine prbnormtndpc(rho,a,b,NDF,N,abseps,IERC,HNC,PRB,BOUND, - * IFAULT) + * IFAULT) double precision A(N),B(N),rho(N),D(N) integer INFIN(N) integer NDF,N,IERC integer IFAULT - double precision HNC,EPS + double precision HNC +C double precision EPS double precision PRB, BOUND double precision, parameter :: infinity = 37.0d0 -Cf2py integer, intent(hide), depend(rho) :: N = len(rho) -Cf2py depend(N) a -Cf2py depend(N) b +Cf2py integer, intent(hide), depend(rho) :: N = len(rho) +Cf2py depend(N) a +Cf2py depend(N) b Cf2py integer, optional :: NDF = 0 Cf2py double precision, optional :: abseps = 0.001 Cf2py double precision, optional :: HNC = 0.24 @@ -30,12 +31,12 @@ CCf2py intent(in) A,B,rho * if INFIN(I) < 0, Ith limits are (-infinity, infinity); * if INFIN(I) = 0, Ith limits are [LOWER(I), infinity); * if INFIN(I) = 1, Ith limits are (-infinity, UPPER(I)]; -* if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)]. +* if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)]. Ndim = 0 DO K = 1,N - Ndim = Ndim + 1 + Ndim = Ndim + 1 INFIN(Ndim) = 2 - D(k) = 0.0 + D(k) = 0.0 if (A(K)-D(K).LE.-INFINITY) THEN if (B(K)-D(K) .GE. INFINITY) THEN Ndim = Ndim - 1 @@ -43,14 +44,14 @@ CCf2py intent(in) A,B,rho else INFIN(Ndim) = 1 endif - else if (B(K)-D(K).GE.INFINITY) THEN + else if (B(K)-D(K).GE.INFINITY) THEN INFIN(Ndim) = 0 endif if (ndim