Directional Spectrum from first 5 moments

[pedrocrxx]

Hi

I’m trying to use the output of the first five moments to recreate the spectra on every valid grid point. I noticed that has been done some minor adjustments to the code (discussion in #103) and I’m still using the latest stable release (v6.07.1), so I’m not sure if my issues are due to that.

When I try to compute a spectrum from the coefficients and compare it against a 2D spectrum directly outputted from the model for the same grid point, I can see some very significant differences between the two versions, with the local peaks on different bins and with different energy values. The overall spectral shape is similar (we can detect the same partitions) but with different peak values (for direction).

I’m doing the following transformation to go from the first five moments to the a1,a2,b1 and b2 coefficients, that later I apply to the MEM:

alpha1=abs(1 - 0.5*(np.radians(sth1m))**2)
alpha2=abs(1 - 2*(np.radians(sth2m))**2)

a1=np.cos(np.radians(th1m))*alpha1
b1=np.sin(np.radians(th1m))*alpha1
a2=np.cos(2*np.radians(th2m))*alpha2
b2=np.sin(2*np.radians(th2m))*alpha2

I’m doing something wrong in this transformation?
I was expecting some differences between the two versions, but not to the extend I’m seeing.

Thank you

Pedro

[mickaelaccensi]

Hi Pedro,

This question is not directly linked to WW3 developer communities, however the computation of the first 5moments to obtain the wave spectrum can be done using the following MEM from Wang He. Notice that your a1,a2, b1, b2 coefficients are correctly computed.

#the MEM method

def buoy_spectrum2d(sf,a1,a2,b1,b2, dirs = np.arange(0,365,10)):
    """(a1,a2,b1,b2,en,ndirs):
    purpose: compute spectrum(f,theta) from coef ai,bi
    author: Wang He
    # Maximum entropy method to estimate the Directional Distribution
    
    # Maximum Entropy Method - Lygre & Krogstad (1986 - JPO)
    # Eqn. 13:
    # phi1 = (c1 - c2c1*)/(1 - abs(c1)^2)
    # phi2 = c2 - c1phi1
    # 2piD = (1 - phi1c1* - phi2c2*)/abs(1 - phi1exp(-itheta) -phi2exp(2itheta))^2
    # c1 and c2 are the complex fourier coefficients
    
    # inumpyuts : dirs in degrees
    """
    nfreq = np.size(sf)
    nbin = np.size(dirs)

    c1 = a1+1j*b1
    c2 = a2+1j*b2
    p1 = (c1-c2*np.conj(c1))/(1.-abs(c1)**2)
    p2 = c2-c1*p1

    # numerator(2D) : x
    x = 1.-p1*np.conj(c1)-p2*np.conj(c2)
    x = numpy.tile(numpy.real(x),(nbin,1)).T

    # denominator(2D): y
    a = dirs*numpy.pi/180.
    e1 = numpy.tile(numpy.cos(a)-1j*numpy.sin(a),(nfreq,1))
    e2 = numpy.tile(numpy.cos(2*a)-1j*numpy.sin(2*a),(nfreq,1))

    p1e1 = numpy.tile(p1,(nbin,1)).T*e1
    p2e2 = numpy.tile(p2,(nbin,1)).T*e2

    y = abs(1-p1e1-p2e2)**2

    D = x/(y*2*numpy.pi)

    # normalizes the spreading function,
    # so that int D(theta,f) dtheta = 1 for each f  

    dth = dirs[1]-dirs[0]

    tot = numpy.tile(numpy.sum(D, axis=1)*dth/180.*numpy.pi,(nbin,1)).T
    D = D/tot

    sp2d = numpy.tile(sf,(nbin,1)).T*D

    return sp2d,D

[pedrocrxx]

I'm applying exactly that same code for the MEM method, so my concern was that I was doing something wrong in the transformation.

I guess that the difference between the directly outputed spectra from ww3_ounp and the spectra computed from the first 5 moments, I posted above, is just the price to pay for the compression.

Thank you for your reply @mickaelaccensi