danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

bugfix in SpecData1D.sim_nl + updated doctests

Browse Source
master
per.andreas.brodtkorb 14 years ago
parent 25906cceb7
commit b93fa2bd68
3 changed files with 88 additions and 60 deletions
Show Stats Download Patch File Download Diff File

110
pywafo/src/wafo/source/c_codes/c_functions.c
Unescape Escape View File

@ -229,64 +229,64 @@ void disufq(double *rvec, double *ivec,
kw1 = kw[ix]; kw1 = kw[ix];
Epij = kw1; Epij = kw1;
for (i=0;i<m;i++,ixi++,iz1++) { for (i=0;i<m;i++,ixi++,iz1++) {
rrA = rA[ixi]*rA[ixi]; /// rrA = rA[ixi]*rA[ixi]; ///
iiA = iA[ixi]*iA[ixi]; /// iiA = iA[ixi]*iA[ixi]; ///
riA = rA[ixi]*iA[ixi]; /// riA = rA[ixi]*iA[ixi]; ///
/// Sum frequency effects along the diagonal /// Sum frequency effects along the diagonal
tmp1 = kfact*(rrA-iiA)*Epij; tmp1 = kfact*(rrA-iiA)*Epij;
tmp2 = kfact*2.0*riA*Epij; tmp2 = kfact*2.0*riA*Epij;
rvec[iz1] += tmp1; rvec[iz1] += tmp1;
ivec[iz1] += tmp2; ivec[iz1] += tmp2;
//rvec[iz2] += tmp1; //rvec[iz2] += tmp1;
//ivec[iz2] -= tmp2; //ivec[iz2] -= tmp2;
//iz2++; //iz2++;
/// Difference frequency effects are zero along the diagonal /// Difference frequency effects are zero along the diagonal
/// and are thus not contributing to the mean. /// and are thus not contributing to the mean.
} }
for (jy = ix+1;jy<nmax;jy++){ for (jy = ix+1;jy<nmax;jy++){
kw2 = kw[jy]; kw2 = kw[jy];
Epij = 0.5*(kw2 + kw1); Epij = 0.5*(kw2 + kw1);
Edij = -0.5*(kw2 - kw1); Edij = -0.5*(kw2 - kw1);
//printf("Edij = %f Epij = %f \n", Edij,Epij); //printf("Edij = %f Epij = %f \n", Edij,Epij);
ixi = ix*m; ixi = ix*m;
jyi = jy*m; jyi = jy*m;
iz1 = ixi+jyi; iz1 = ixi+jyi;
iv1 = jyi-ixi; iv1 = jyi-ixi;
//iz2 = (n*m-iz1); //iz2 = (n*m-iz1);
//iv2 = (n*m-iv1); //iv2 = (n*m-iv1);
for (i = 0;i<m;i++,ixi++,jyi++,iz1++,iv1++) { for (i = 0;i<m;i++,ixi++,jyi++,iz1++,iv1++) {
rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy]; rrA = rA[ixi]*rA[jyi]; ///rrA = rA[i][ix]*rA[i][jy];
iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy]; iiA = iA[ixi]*iA[jyi]; ///iiA = iA[i][ix]*iA[i][jy];
riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy]; riA = rA[ixi]*iA[jyi]; ///riA = rA[i][ix]*iA[i][jy];
irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy]; irA = iA[ixi]*rA[jyi]; ///irA = iA[i][ix]*rA[i][jy];
/* Sum frequency effects */ /* Sum frequency effects */
tmp1 = kfact*2.0*(rrA-iiA)*Epij; tmp1 = kfact*2.0*(rrA-iiA)*Epij;
tmp2 = kfact*2.0*(riA+irA)*Epij; tmp2 = kfact*2.0*(riA+irA)*Epij;
rvec[iz1] += tmp1;///rvec[i][ix+jy] += tmp1; rvec[iz1] += tmp1;///rvec[i][ix+jy] += tmp1;
ivec[iz1] += tmp2;///ivec[i][ix+jy] += tmp2; ivec[iz1] += tmp2;///ivec[i][ix+jy] += tmp2;
//rvec[iz2] += tmp1;///rvec[i][n*m-(ix+jy)] += tmp1; //rvec[iz2] += tmp1;///rvec[i][n*m-(ix+jy)] += tmp1;
//ivec[iz2] -= tmp2;///ivec[i][n*m-(ix+jy)] -= tmp2; //ivec[iz2] -= tmp2;///ivec[i][n*m-(ix+jy)] -= tmp2;
// iz2++; // iz2++;
/* Difference frequency effects */ /* Difference frequency effects */
tmp1 = kfact*2.0*(rrA+iiA)*Edij; tmp1 = kfact*2.0*(rrA+iiA)*Edij;
tmp2 = kfact*2.0*(riA-irA)*Edij; tmp2 = kfact*2.0*(riA-irA)*Edij;
rvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1; rvec[iv1] += tmp1;///rvec[i][jy-ix] += tmp1;
ivec[iv1] += tmp2;///ivec[i][jy-ix] += tmp2; ivec[iv1] += tmp2;///ivec[i][jy-ix] += tmp2;
//rvec[iv2] += tmp1;///rvec[i][n*m-(jy-ix)] += tmp1; //rvec[iv2] += tmp1;///rvec[i][n*m-(jy-ix)] += tmp1;
//ivec[iv2] -= tmp2;///ivec[i][n*m-(jy-ix)] -= tmp2; //ivec[iv2] -= tmp2;///ivec[i][n*m-(jy-ix)] -= tmp2;
//iv2++; //iv2++;
} }
} }
} }
} }
else{ /* Finite water depth */ else{ /* Finite water depth */
for (ix = nmin-1;ix<nmax;ix++) { for (ix = nmin-1;ix<nmax;ix++) {

50
pywafo/src/wafo/spectrum/core.py
Unescape Escape View File

@ -1496,15 +1496,11 @@ class SpecData1D(WafoData):
... res = fun(x2[:,1::],axis=0) ... res = fun(x2[:,1::],axis=0)
... m = res.mean() ... m = res.mean()
... sa = res.std() ... sa = res.std()
... trueval, m, sa ... #trueval, m, sa
... np.abs(m-trueval)<sa ... np.abs(m-trueval)<sa
(0, 1.2309597785531439e-19, 3.2562166904166163e-18)
True True
(array([ 1.74952003]), 1.7500502911292997, 0.022461280887884415)
array([ True], dtype=bool) array([ True], dtype=bool)
(0.0, -0.00024040615507690482, 0.0087615749174770451)
True True
(0.0, 0.0018609047569154713, 0.049873521257196997)
True True
waveplot(x1,'r',x2,'g',1,1) waveplot(x1,'r',x2,'g',1,1)
@ -1726,7 +1722,29 @@ class SpecData1D(WafoData):
... res = fun(x2[:,1::], axis=0) ... res = fun(x2[:,1::], axis=0)
... m = res.mean() ... m = res.mean()
... sa = res.std() ... sa = res.std()
... trueval, m, sa ... #trueval, m, sa
... np.abs(m-trueval)<sa
True
True
True
True
>>> x = []
>>> for i in range(20):
... x2, x1 = S.sim_nl(ns=20000,cases=1)
... x.append(x2[:,1::])
>>> x2 = np.hstack(x)
>>> truth1 = [0,np.sqrt(S.moment(1)[0][0])] + S.stats_nl(moments='sk')
>>> truth1[-1] = truth1[-1]-3
>>> truth1
[0, 1.7495200310090628, 0.18673120577479821, 0.06198852126241805]
>>> funs = [np.mean,np.std,st.skew,st.kurtosis]
>>> for fun,trueval in zip(funs,truth1):
... res = fun(x2, axis=0)
... m = res.mean()
... sa = res.std()
... #trueval, m, sa
... np.abs(m-trueval)<sa ... np.abs(m-trueval)<sa
True True
True True
@ -1889,10 +1907,20 @@ class SpecData1D(WafoData):
## % 1'st order + 2'nd order component. ## % 1'st order + 2'nd order component.
## x2(:,2:end) =x(:,2:end)+ real(x2s(1:np,:))+real(x2d(1:np,:)) ## x2(:,2:end) =x(:,2:end)+ real(x2s(1:np,:))+real(x2d(1:np,:))
## else ## else
amp = np.array(amp.T).ravel() if False:
rvec, ivec = c_library.disufq(amp.real, amp.imag, w, kw, water_depth, g, nmin, nmax, cases, ns) # TODO: disufq does not work for cases>1
amp = np.array(amp.T).ravel()
svec = rvec + 1J * ivec rvec, ivec = c_library.disufq(amp.real, amp.imag, w, kw, water_depth,
g, nmin, nmax, cases, ns)
svec = rvec + 1J * ivec
else:
amp = amp.T
svec=[]
for i in range(cases):
rvec, ivec = c_library.disufq(amp[i].real, amp[i].imag, w, kw, water_depth,
g, nmin, nmax, 1, ns)
svec.append(rvec + 1J * ivec)
svec = np.hstack(svec)
svec.shape = (cases, ns) svec.shape = (cases, ns)
x2o = fft(svec, axis=1).T # 2'nd order component x2o = fft(svec, axis=1).T # 2'nd order component
@ -3022,7 +3050,7 @@ class SpecData2D(WafoData):
Example: Example:
>>> import wafo.spectrum.models as sm >>> import wafo.spectrum.models as sm
>>> D = sm.Spreading() >>> D = sm.Spreading()
>>> SD = D.tospecdata2d(sm.Jonswap().tospecdata()) >>> SD = D.tospecdata2d(sm.Jonswap().tospecdata(),nt=101)
>>> m,mtext = SD.moment(nr=2,vari='xyt') >>> m,mtext = SD.moment(nr=2,vari='xyt')
>>> np.round(m,3),mtext >>> np.round(m,3),mtext
(array([ 3.061, 0.132, -0. , 2.13 , 0.011, 0.008, 1.677, -0. , (array([ 3.061, 0.132, -0. , 2.13 , 0.011, 0.008, 1.677, -0. ,

6
pywafo/src/wafo/spectrum/models.py
Unescape Escape View File

@ -447,9 +447,9 @@ class Jonswap(ModelSpectrum):
N : scalar defining decay of high frequency part. (default 5) N : scalar defining decay of high frequency part. (default 5)
M : scalar defining spectral width around the peak. (default 4) M : scalar defining spectral width around the peak. (default 4)
method : String defining method used to estimate Ag when gamma>1 method : String defining method used to estimate Ag when gamma>1
'integrate' : Ag = 1/gaussq(Gf*ggamspec(wn,N,M),0,wnc) (default) 'integration': Ag = 1/gaussq(Gf*ggamspec(wn,N,M),0,wnc) (default)
'parametric': Ag = (1+f1(N,M)*log(gamma)**f2(N,M))/gamma 'parametric' : Ag = (1+f1(N,M)*log(gamma)**f2(N,M))/gamma
'custom' : Ag = Ag 'custom' : Ag = Ag
wnc : wc/wp normalized cut off frequency used when calculating Ag wnc : wc/wp normalized cut off frequency used when calculating Ag
by integration (default 6) by integration (default 6)
Parameters Parameters

Write Preview
Loading…
Cancel
Save