From 2e3e506767a5c220223e91a48dc30390fcebd26a Mon Sep 17 00:00:00 2001 From: Per A Brodtkorb Date: Fri, 16 Mar 2018 19:53:57 +0100 Subject: [PATCH] Added extra tests for dispersion_relation.py and qtf in spectrum/core.py. Deleted powerpoint.py --- wafo/objects.py | 18 +- wafo/powerpoint.py | 311 ------------------ wafo/spectrum/core.py | 217 +++++++++--- wafo/spectrum/models.py | 2 +- wafo/wave_theory/dispersion_relation.py | 38 +-- .../tests/test_dispersion_relation.py | 45 ++- 6 files changed, 245 insertions(+), 386 deletions(-) delete mode 100644 wafo/powerpoint.py diff --git a/wafo/objects.py b/wafo/objects.py index cd04f5b..0d09ac8 100644 --- a/wafo/objects.py +++ b/wafo/objects.py @@ -762,6 +762,7 @@ class CyclePairs(PlotData): >>> mm = tp.cycle_pairs() >>> lc = mm.level_crossings() + h = mm.plot(marker='.') h2 = lc.plot() @@ -959,8 +960,17 @@ class CyclePairs(PlotData): F = Estimated cycle matrix. h = Selected bandwidth. - See also - dcc2cmat, cc2dcc, smoothcmat + Example: + -------- + >>> import wafo + >>> ts = wafo.objects.mat2timeseries(wafo.data.sea()) + >>> tp = ts.turning_points() + >>> mm = tp.cycle_pairs() + >>> cm = mm.cycle_matrix((-3, 3, 50)) + + + See also + dcc2cmat, cc2dcc, smoothcmat """ if not 0 <= method <= 2: @@ -980,7 +990,7 @@ class CyclePairs(PlotData): # if method >= 1: # F, h = smoothcmat(F,method, h, NOsubzero, alpha) - return CycleMatrix(F, u, u) + return CycleMatrix(F, (u, u)) def _dcp2cmat(self, dcp, n): """ @@ -1065,7 +1075,7 @@ class CyclePairs(PlotData): if np.any(ix): cp1[ix], cp2[ix] = cp2[ix], cp1[ix] - return np.asarray(cp1, type=int), np.asarray(cp2, type=int) + return np.asarray(cp1, dtype=int), np.asarray(cp2, dtype=int) class TurningPoints(PlotData): diff --git a/wafo/powerpoint.py b/wafo/powerpoint.py deleted file mode 100644 index ffa151a..0000000 --- a/wafo/powerpoint.py +++ /dev/null @@ -1,311 +0,0 @@ -''' -Created on 15. des. 2009 - -@author: pab -''' -# import os -# import sys -# import win32com -# from win32com.client.selecttlb import EnumTlbs -# typelib_mso = None -# typelib_msppt = None -# for typelib in EnumTlbs(): -# d = typelib.desc.split(' ') -# if d[0] == 'Microsoft' and d[1] == 'Office' and d[3] == 'Object' \ -# and d[4] == 'Library': -# typelib_mso = typelib -# if d[0] == 'Microsoft' and d[1] == 'PowerPoint' and d[3] == 'Object' \ -# and d[4] == 'Library': -# typelib_msppt = typelib -# if hasattr(sys, 'frozen'): # If we're an .exe file -# win32com.__gen_path__ = os.path.dirname(sys.executable) -# # win32com.__gen_path__ = os.environ['TEMP'] -# if win32com.client.gencache.is_readonly: -# win32com.client.gencache.is_readonly = False -# win32com.client.gencache.Rebuild() -# MSPPT = win32com.client.gencache.EnsureModule(typelib_msppt.clsid, -# typelib_msppt.lcid, -# int(typelib_msppt.major), -# int(typelib_msppt.minor)) -# MSO = win32com.client.gencache.EnsureModule(typelib_mso.clsid, -# typelib_mso.lcid, -# int(typelib_mso.major), int(typelib_mso.minor)) -from __future__ import absolute_import -from six import iteritems # @UnresolvedImport -import os -import warnings -import win32com.client -from . import MSO -from . import MSPPT -from PIL import Image # @UnresolvedImport - -g = globals() -for c in dir(MSO.constants): - g[c] = getattr(MSO.constants, c) -for c in dir(MSPPT.constants): - g[c] = getattr(MSPPT.constants, c) - - -class Powerpoint(object): - - def __init__(self, file_name=''): - - self.application = win32com.client.Dispatch("Powerpoint.Application") - # self.application.Visible = True - self._visible = self.application.Visible - if file_name: - self.presentation = self.application.Presentations.Open(file_name) - else: - self.presentation = self.application.Presentations.Add() - self.num_slides = 0 - # default picture width and height - self.default_width = 500 - self.default_height = 400 - self.title_font = 'Arial' # 'Boopee' - self.title_size = 36 - self.text_font = 'Arial' # 'Boopee' - self.text_size = 20 - self.footer = '' - - def set_footer(self): - ''' - Set Footer in SlideMaster and NotesMaster - ''' - if self.footer: - if self.presentation.HasTitleMaster: - TMHF = self.presentation.TitleMaster.HeadersFooters - TMHF.Footer.Text = self.footer - TMHF.Footer.Visible = True - - SMHF = self.presentation.SlideMaster.HeadersFooters - SMHF.Footer.Text = self.footer - SMHF.Footer.Visible = True - SMHF.SlideNumber.Visible = True - NMHF = self.presentation.NotesMaster.HeadersFooters - NMHF.Footer.Text = self.footer - NMHF.SlideNumber.Visible = True - for slide in self.presentation.Slides: - shapes = slide.Shapes - for shape in shapes: - if shape.Name == 'Footer': - footer = shape - break - else: - footer = shapes.AddTextbox( - msoTextOrientationHorizontal, # @UndefinedVariable - Left=0, Top=510, Width=720, Height=28.875) - footer.Name = 'Footer' - footer.TextFrame.TextRange.Text = self.footer - - def add_title_slide(self, title, subtitle=''): - self.num_slides += 1 - slide = self.presentation.Slides.Add( - self.num_slides, MSPPT.constants.ppLayoutTitle) - - unused_title_id, unused_textbox_id = 1, 2 - for id_, title1 in enumerate([title, subtitle]): - titlerange = slide.Shapes(id_ + 1).TextFrame.TextRange - titlerange.Text = title1 - titlerange.Font.Name = self.title_font - titlerange.Font.Size = self.title_size - id_ * \ - 12 if self.title_size > 22 else self.title_size - - def add_slide(self, title='', texts='', notes='', image_file='', - maxlevel=None, left=220, width=-1, height=-1): - self.num_slides += 1 - slide = self.presentation.Slides.Add( - self.num_slides, MSPPT.constants.ppLayoutText) - - self.add2slide(slide, title, texts, notes, image_file, maxlevel, left, - width, height) - return slide - - def add2slide(self, slide, title='', texts='', notes='', image_file='', - maxlevel=None, left=220, width=-1, height=-1, - keep_aspect=True): - title_id, textbox_id = 1, 2 - if title: - titlerange = slide.Shapes(title_id).TextFrame.TextRange - titlerange.Font.Name = self.title_font - titlerange.Text = title - titlerange.Font.Size = self.title_size - - if texts != '' and texts != ['']: - # textrange = slide.Shapes(textbox_id).TextFrame.TextRange - self._add_text(slide, textbox_id, texts, maxlevel) - - if image_file != '' and image_file != ['']: - if keep_aspect: - im = Image.open(image_file) - t_w, t_h = im.size - if height <= 0 and width <= 0: - if t_w * self.default_height < t_h * self.default_width: - height = self.default_height - else: - width = self.default_width - if height <= 0 and width: - height = t_h * width / t_w - elif height and width <= 0: - width = t_w * height / t_h - - slide.Shapes.AddPicture(FileName=image_file, LinkToFile=False, - SaveWithDocument=True, - Left=left, Top=110, - Width=width, Height=height) # 400) - if notes != '' and notes != ['']: - notespage = slide.NotesPage # .Shapes(2).TextFrame.TextRange - self._add_text(notespage, 2, notes) - return slide - - def _add_text(self, page, id, txt, maxlevel=None): # @ReservedAssignment - page.Shapes(id).TextFrame.TextRange.Font.Name = self.text_font - - if isinstance(txt, dict): - self._add_text_from_dict(page, id, txt, 1, maxlevel) - elif isinstance(txt, (list, tuple)): - self._add_text_from_list(page, id, txt, maxlevel) - else: - unused_tr = page.Shapes(id).TextFrame.TextRange.InsertAfter(txt) - unused_temp = page.Shapes(id).TextFrame.TextRange.InsertAfter('\r') - - page.Shapes(id).TextFrame.TextRange.Font.Size = self.text_size - - def _add_text_from_dict(self, page, id, txt_dict, # @ReservedAssignment - level, maxlevel=None): - if maxlevel is None or level <= maxlevel: - for name, subdict in iteritems(txt_dict): - tr = page.Shapes(id).TextFrame.TextRange.InsertAfter(name) - unused_temp = page.Shapes( - id).TextFrame.TextRange.InsertAfter('\r') - tr.IndentLevel = level - self._add_text_from_dict( - page, id, subdict, min(level + 1, 5), maxlevel) - - def _add_text_from_list(self, page, id, # @ReservedAssignment - txt_list, maxlevel=None): - for txt in txt_list: - level = 1 - while isinstance(txt, (list, tuple)): - txt = txt[0] - level += 1 - if maxlevel is None or level <= maxlevel: - tr = page.Shapes(id).TextFrame.TextRange.InsertAfter(txt) - unused_temp = page.Shapes( - id).TextFrame.TextRange.InsertAfter('\r') - tr.IndentLevel = level - - def save(self, fullfile=''): - if fullfile: - self.presentation.SaveAs(FileName=fullfile) - else: - self.presentation.Save() - - def quit(self): # @ReservedAssignment - if self._visible: - self.presentation.Close() - else: - self.application.Quit() - - def quit_only_if_hidden(self): - if not self._visible: - self.application.Quit() - - -def test_powerpoint(): - # Make powerpoint - - ppt = Powerpoint() - ppt.footer = 'This is the footer' - ppt.add_title_slide('Title', 'Per A.') - ppt.add_slide(title='alsfkasldk', texts='asdflaf', notes='asdfas') - ppt.set_footer() - - -def make_ppt(): - application = win32com.client.Dispatch("Powerpoint.Application") - application.Visible = True - presentation = application.Presentations.Add() - slide1 = presentation.Slides.Add(1, MSPPT.constants.ppLayoutText) - -# title = slide1.Shapes.AddTextBox(Type=msoTextOrientationHorizontal, -# Left=50, Top=10, Width=620, Height=70) -# title.TextFrame.TextRange.Text = 'Overskrift' - title_id, textbox_id = 1, 2 - slide1.Shapes(title_id).TextFrame.TextRange.Text = 'Overskrift' - # slide1.Shapes(title_id).TextFrame.Width = 190 - - slide1.Shapes(textbox_id).TextFrame.TextRange.InsertAfter('Test') - unused_tr = slide1.Shapes(textbox_id).TextFrame.TextRange.InsertAfter('\r') - slide1.Shapes(textbox_id).TextFrame.TextRange.IndentLevel = 1 - tr = slide1.Shapes(textbox_id).TextFrame.TextRange.InsertAfter('tests') - unused_tr0 = slide1.Shapes( - textbox_id).TextFrame.TextRange.InsertAfter('\r') - tr.IndentLevel = 2 - tr1 = slide1.Shapes(textbox_id).TextFrame.TextRange.InsertAfter('test3') - tr1.IndentLevel = 3 - # slide1.Shapes(textbox_id).TextFrame.TextRange.Text = 'Test \r test2' - -# textbox = slide1.Shapes.AddTextBox(Type=msoTextOrientationHorizontal, -# Left=30, Top=100, Width=190, Height=400) -# textbox.TextFrame.TextRange.Text = 'Test \r test2' - # picbox = slide1.Shapes(picb_id) - - filename = r'c:\temp\data1_report1_and_2_Tr120_1.png' - slide1.Shapes.AddPicture(FileName=filename, LinkToFile=False, - SaveWithDocument=True, - Left=220, Top=100, Width=500, Height=420) - - slide1.NotesPage.Shapes(2).TextFrame.TextRange.Text = 'test' - - -# for shape in slide1.Shapes: -# shape.TextFrame.TextRange.Text = 'Test \r test2' - # slide1.Shapes.Titles.TextFrames.TestRange.Text -# shape = slide1.Shapes.AddShape(msoShapeRectangle, 300, 100, 400, 400) -# shape.TextFrame.TextRange.Text = 'Test \n test2' -# shape.TextFrame.TextRange.Font.Size = 12 - # -# app = wx.PySimpleApp() -# dialog = wx.FileDialog(None, 'Choose image file', defaultDir=os.getcwd(), -# wildcard='*.*', -# style=wx.OPEN | wx.CHANGE_DIR | wx.MULTIPLE) -# -# if dialog.ShowModal() == wx.ID_OK: -# files_or_paths = dialog.GetPaths() -# for filename in files_or_paths: -# slide1.Shapes.AddPicture(FileName=filename, LinkToFile=False, -# SaveWithDocument=True, -# Left=100, Top=100, Width=200, Height=200) -# dialog.Destroy() - # presentation.Save() - # application.Quit() -def rename_ppt(): - root = r'C:/pab/tsm_opeval/analysis_tsmps_aco_v2008b/plots' - filenames = os.listdir(root) - prefix = 'TSMPSv2008b_' - for filename in filenames: - if filename.endswith('.ppt'): - try: - ppt = Powerpoint(os.path.join(root, filename)) - ppt.footer = prefix + filename - ppt.set_footer() - ppt.save(os.path.join(root, ppt.footer)) - except: - warnings.warn('Unable to load %s' % filename) - - -def load_file_into_ppt(): - root = r'C:/pab/tsm_opeval/analysis_tsmps_aco_v2008b/plots' - filenames = os.listdir(root) - prefix = 'TSMPSv2008b_' - for filename in filenames: - if filename.startswith(prefix) and filename.endswith('.ppt'): - try: - unused_ppt = Powerpoint(os.path.join(root, filename)) - except: - warnings.warn('Unable to load %s' % filename) -if __name__ == '__main__': - # make_ppt() - # test_powerpoint() - # load_file_into_ppt() - rename_ppt() diff --git a/wafo/spectrum/core.py b/wafo/spectrum/core.py index b4b0590..402f819 100644 --- a/wafo/spectrum/core.py +++ b/wafo/spectrum/core.py @@ -50,6 +50,7 @@ __all__ = ['SpecData1D', 'SpecData2D', 'plotspec'] _EPS = np.finfo(float).eps +_TINY = np.finfo(float).tiny def _set_seed(iseed): @@ -61,7 +62,7 @@ def _set_seed(iseed): random.seed(iseed) -def qtf(w, h=inf, g=9.81): +def qtf(w, h=inf, g=9.81, method='winterstein', rtol=1e-7, atol=0): """ Return Quadratic Transfer Function @@ -79,69 +80,212 @@ def qtf(w, h=inf, g=9.81): h_s = sum frequency effects h_d = difference frequency effects h_dii = diagonal of h_d + + Example + ------- + + >>> w = np.r_[0.1, 1./3, 2./3, 1] + >>> hs, hd, hdi = qtf(w, h=np.inf, g=9.81) + >>> np.allclose(hs, [[ 0.00050968, 0.00308642, 0.01158115, 0.02573904], + ... [ 0.00308642, 0.00566316, 0.01415789, 0.02831578], + ... [ 0.01158115, 0.01415789, 0.02265262, 0.03681051], + ... [ 0.02573904, 0.02831578, 0.03681051, 0.0509684 ]]) + True + + >>> np.allclose(hd, [[-0. , -0.00257674, -0.01107147, -0.02522936], + ... [-0.00257674, -0. , -0.00849473, -0.02265262], + ... [-0.01107147, -0.00849473, -0. , -0.01415789], + ... [-0.02522936, -0.02265262, -0.01415789, -0. ]]) + True + + >>> hs2, hd2, hdi2 = qtf(w, h=1e+6, g=9.81, method='winterstein') + >>> np.allclose(hs2, [[0.00050968, 0.00308642, 0.01158115, 0.02573904], + ... [0.00308642, 0.00566316, 0.01415789, 0.02831578], + ... [0.01158115, 0.01415789, 0.02265262, 0.03681051], + ... [0.02573904, 0.02831578, 0.03681051, 0.0509684 ]]) + True + >>> np.allclose(hd2, [[-2.50061328e-07, 1.38729557e-03, 8.18314621e-03, 2.06421189e-02], + ... [1.38729557e-03, -2.50005518e-07, 2.83135545e-03, 1.13261230e-02], + ... [8.18314621e-03, 2.83135545e-03, -2.50001380e-07, 2.83133750e-03], + ... [2.06421189e-02, 1.13261230e-02, 2.83133750e-03, -2.50000613e-07]]) + True + + >>> w = np.r_[0, 1e-10, 1e-5, 1e-1] + >>> hs, hd, hdi = qtf(w, h=np.inf, g=9.81) + >>> np.allclose(hs, [[0.00000000e+00, 2.54841998e-22, 2.54841998e-12, 2.54841998e-04], + ... [2.54841998e-22, 5.09683996e-22, 2.54841998e-12, 2.54841998e-04], + ... [2.54841998e-12, 2.54841998e-12, 5.09683996e-12, 2.54842001e-04], + ... [2.54841998e-04, 2.54841998e-04, 2.54842001e-04, 5.09683996e-04]]) + True + + >>> np.allclose(hd, [[-0.00000000e+00, -2.54841998e-22, -2.54841998e-12, -2.54841998e-04], + ... [-2.54841998e-22, -0.00000000e+00, -2.54841998e-12, -2.54841998e-04], + ... [-2.54841998e-12, -2.54841998e-12, -0.00000000e+00, -2.54841995e-04], + ... [-2.54841998e-04, -2.54841998e-04, -2.54841995e-04, -0.00000000e+00]]) + True + + >>> hs2, hd2, hdi2 = qtf(w, h=1e+100, g=9.81, method='winterstein') + >>> np.allclose(hs2, [[1.50234572e-63, 2.54841998e-22, 2.54841998e-12, 2.54841998e-04], + ... [2.54841998e-22, 5.09683996e-22, 2.54841998e-12, 2.54841998e-04], + ... [2.54841998e-12, 2.54841998e-12, 5.09683996e-12, 2.54842001e-04], + ... [2.54841998e-04, 2.54841998e-04, 2.54842001e-04, 5.09683996e-04]], + ... atol=0) + True + + >>> np.allclose(hd2, [[-2.50000000e-101, 2.54841998e-022, 2.54841998e-012, 2.54841998e-004], + ... [2.54841998e-022, -2.50000000e-101, 2.54836901e-012, 2.54841997e-004], + ... [2.54841998e-012, 2.54836901e-012, -2.50000000e-101, 2.54791032e-004], + ... [2.54841998e-004, 2.54841997e-004, 2.54791032e-004, -2.50000000e-101]], + ... atol=0) + True + + References: + ----------- + Langley, RS (1987) + 'A statistical analysis of nonlinear random waves' + Ocean Engineering, Vol 14, No 5, pp 389-407 + + Marthinsen, T. and Winterstein, S.R (1992) + 'On the skewness of random surface waves' + In proceedings of the 2nd ISOPE Conference, San Francisco, 14-19 june. """ +# >>> hs3, hd3, hdi3 = qtf(w, h=200, g=9.81, method='winterstein') +# >>> hs3 +# +# >>> hd3 +# +# >>> np.allclose(hs3, [[ 0. , 0.00283158, 0.01132631, 0.0254842 ], +# ... [ 0.00283158, 0.00566316, 0.01415789, 0.02831578], +# ... [ 0.01132631, 0.01415789, 0.02265262, 0.03681051], +# ... [ 0.0254842 , 0.02831578, 0.03681051, 0.0509684 ]]) +# +# >>> np.allclose(hd3, [[-0. , -0.00283158, -0.01132631, -0.0254842 ], +# ... [-0.00283158, -0. , -0.00849473, -0.02265262], +# ... [-0.01132631, -0.00849473, -0. , -0.01415789], +# ... [-0.0254842 , -0.02265262, -0.01415789, -0. ]]) + w = atleast_1d(w) num_w = w.size - k_w = w2k(w, theta=0, h=h, g=g)[0] - - k_1, k_2 = meshgrid(k_w, k_w) - if h == inf: # go here for faster calculations + k_w = w2k(w, theta=0, h=h, g=g, rtol=rtol, atol=atol)[0] + k_1, k_2 = meshgrid(k_w, k_w, sparse=True) h_s = 0.25 * (abs(k_1) + abs(k_2)) h_d = -0.25 * abs(abs(k_1) - abs(k_2)) h_dii = zeros(num_w) return h_s, h_d, h_dii - [w_1, w_2] = meshgrid(w, w) - - w12 = (w_1 * w_2) - w1p2 = (w_1 + w_2) - w1m2 = (w_1 - w_2) - k12 = (k_1 * k_2) - k1p2 = (k_1 + k_2) + w1 = w + _TINY **(1./10) * (np.sign(w) * np.int_(np.abs(w) < _EPS) + np.int_(w==0)) + # dw = w1-w + # print(dw) + w = w1 + # k_w += _TINY ** (1./3) * (np.sign(k_w) * np.int_(np.abs(k_w) < _EPS) + np.int_(k_w==0)) + k_w = w2k(w, theta=0, h=h, g=g, rtol=rtol, atol=atol)[0] + k_1, k_2 = meshgrid(k_w, k_w, sparse=True) + w_1, w_2 = meshgrid(w, w, sparse=True) + + w12 = w_1 * w_2 + w1p2 = w_1 + w_2 + w1m2 = w_1 - w_2 + k12 = k_1 * k_2 + k1p2 = k_1 + k_2 k1m2 = abs(k_1 - k_2) - if False: # Langley + if method.startswith('langley'): p_1 = (-2 * w1p2 * (k12 * g ** 2. - w12 ** 2.) + w_1 * (w_2 ** 4. - g ** 2 * k_2 ** 2) + - w_2 * (w_1 ** 4 - g * 2. * k_1 ** 2)) / (4. * w12) + w_2 * (w_1 ** 4 - g * 2. * k_1 ** 2)) / (4. * w12 +_TINY) + p_2 = w1p2 ** 2. * cosh((k1p2) * h) - g * (k1p2) * sinh((k1p2) * h) + p_2 += _TINY * np.int_(p_2==0) h_s = (-p_1 / p_2 * w1p2 * cosh((k1p2) * h) / g - - (k12 * g ** 2 - w12 ** 2.) / (4 * g * w12) + + (k12 * g ** 2 - w12 ** 2.) / (4 * g * w12 + _TINY) + (w_1 ** 2 + w_2 ** 2) / (4 * g)) p_3 = (-2 * w1m2 * (k12 * g ** 2 + w12 ** 2) - w_1 * (w_2 ** 4 - g ** 2 * k_2 ** 2) + - w_2 * (w_1 ** 4 - g ** 2 * k_1 ** 2)) / (4. * w12) + w_2 * (w_1 ** 4 - g ** 2 * k_1 ** 2)) / (4. * w12 + _TINY) p_4 = w1m2 ** 2. * cosh(k1m2 * h) - g * (k1m2) * sinh((k1m2) * h) + p_4 += _TINY * np.int_(p_4==0) h_d = (-p_3 / p_4 * (w1m2) * cosh((k1m2) * h) / g - - (k12 * g ** 2 + w12 ** 2) / (4 * g * w12) + + (k12 * g ** 2 + w12 ** 2) / (4 * g * w12 + _TINY) + (w_1 ** 2. + w_2 ** 2.) / (4. * g)) else: # Marthinsen & Winterstein - tmp1 = 0.5 * g * k12 / w12 - tmp2 = 0.25 / g * (w_1 ** 2. + w_2 ** 2. + w12) - h_s = (tmp1 - tmp2 + 0.25 * g * (w_1 * k_2 ** 2. + w_2 * k_1 ** 2) / - (w12 * (w1p2))) / (1. - g * (k1p2) / (w1p2) ** 2. * - tanh((k1p2) * h)) + tmp2 - 0.5 * tmp1 # OK - - tmp2 = 0.25 / g * (w_1 ** 2 + w_2 ** 2 - w12) # OK - h_d = (tmp1 - tmp2 - 0.25 * g * (w_1 * k_2 ** 2 - w_2 * k_1 ** 2) / - (w12 * (w1m2))) / (1. - g * (k1m2) / (w1m2) ** 2. * - tanh((k1m2) * h)) + tmp2 - 0.5 * tmp1 # OK - - # tmp1 = 0.5*g*k_w./(w.*sqrt(g*h)) - # tmp2 = 0.25*w.^2/g + tmp1 = 2.0 * g * k12 / (w12 + 0) + tmp2 = (w_1 ** 2. + w_2 ** 2. + w12) / g + + h_s = 0.25 * ((tmp1 - tmp2 + + g * (w_1 * k_2 ** 2. + w_2 * k_1 ** 2) / (w12 * w1p2 + 0)) + / (1. - g * h * k1p2 / (w1p2 ** 2. + 0) * tanh(k1p2)) + + tmp2 - 0.5 * tmp1) # OK + + tiny_diag = _TINY * np.diag(np.ones(num_w)) # Avoid division by zero on diagonal + tmp3 = (w_1 ** 2 + w_2 ** 2 - w12) / g # OK + numerator = (tmp1 - tmp3 - g * (w_1 * k_2 ** 2 - w_2 * k_1 ** 2) / (w12 * w1m2 + tiny_diag)) + h_d = 0.25 * (numerator / (1. - g * h * k1m2 / (w1m2 ** 2. + tiny_diag) * tanh(k1m2)) + + tmp3 - 0.5 * tmp1) # OK + +# h_d = 0.25 * ((tmp1 - tmp3 +# - g * (w_1 * k_2 ** 2 - w_2 * k_1 ** 2) / (w12 * w1m2 + tiny_diag)) +# / (1. - g * h * k1m2 / (w1m2 ** 2. + tiny_diag) * tanh(k1m2)) +# + tmp3 - 0.5 * tmp1) # OK + + # tmp1 = 2 * g * (k_w./w)^2 + # tmp2 = w.^2/g # Wave group velocity - c_g = 0.5 * g * (tanh(k_w * h) + k_w * h * (1.0 - tanh(k_w * h) ** 2)) / w - h_dii = (0.5 * (0.5 * g * (k_w / w) ** 2. - 0.5 * w ** 2 / g + - g * k_w / (w * c_g)) / - (1. - g * h / c_g ** 2.) - 0.5 * k_w / sinh(2 * k_w * h)) # OK + k_h = k_w * h + c_g = 0.5 * g * (tanh(k_h) + k_h / np.cosh(k_h) ** 2) / w + numerator2 = (g * (k_w / w) ** 2. - w ** 2 / g + 2 * g * k_w / (w * c_g + 0)) + h_dii = 0.25 * (numerator2 / (1. - g * h / (c_g ** 2. + 0)) + - 2 * k_w / sinh(2 * k_h)) # OK +# c_g = 0.5 * g * (tanh(k_w * h) + k_w * h * (1.0 - tanh(k_w * h) ** 2)) / w +# h_dii = (0.5 * (0.5 * g * (k_w / w) ** 2. - 0.5 * w ** 2 / g + g * k_w / (w * c_g + 0)) +# / (1. - g * h / (c_g ** 2. + 0)) +# - 0.5 * k_w / sinh(2 * k_w * h)) # OK h_d.flat[0::num_w + 1] = h_dii + # infinite water + # >>> np.allclose(hs, [[ 0. , 0.00283158, 0.01132631, 0.0254842 ], + # ... [ 0.00283158, 0.00566316, 0.01415789, 0.02831578], + # ... [ 0.01132631, 0.01415789, 0.02265262, 0.03681051], + # ... [ 0.0254842 , 0.02831578, 0.03681051, 0.0509684 ]]) + # True + # + # >>> np.allclose(hd, [[-0. , -0.00283158, -0.01132631, -0.0254842 ], + # ... [-0.00283158, -0. , -0.00849473, -0.02265262], + # ... [-0.01132631, -0.00849473, -0. , -0.01415789], + # ... [-0.0254842 , -0.02265262, -0.01415789, -0. ]]) + # True + + # winterstein + # h_s = + # [[ 0.00000000e+00 -1.64775418e+00 -6.95612056e-01 -4.18817231e-01 -3.15690232e-01] + # [ -1.64775418e+00 -7.98574421e-01 -2.21051428e-01 -6.73808482e-02 -1.69373060e-02] + # [ -6.95612056e-01 -2.21051428e-01 -1.29139936e-01 -4.11797418e-02 1.12063541e-03] + # [ -4.18817231e-01 -6.73808482e-02 -4.11797418e-02 -3.51718594e-03 2.44489725e-02] + + # h_d = + # [[ 0. -1.64775418 -0.69561206 -0.41881723 -0.31569023] + # [-1.6103978 0.0130494 -0.12128861 -0.05785645 -0.02806048] + # [-0.65467824 -0.12128861 0.01494093 -0.04402996 -0.02595442] + # [-0.35876732 -0.05785645 -0.04402996 0.01905565 -0.02218373] + # [ inf -0.02806048 -0.02595442 -0.02218373 0.02705736]] + # langley + # h_d = [[ 0.00000000e+00 -8.87390092e+14 -3.87924869e+14 -1.66844106e+15 inf] + # [ -8.87390092e+14 -1.14566397e-02 -1.50113192e-01 -1.11791139e-01 -1.13090565e-01] + # [ -3.87924869e+14 -1.50113192e-01 -8.56987798e-03 -5.10233013e-02 -4.93936523e-02] + # [ -1.66844106e+15 -1.11791139e-01 -5.10233013e-02 -4.72078473e-03 -2.74040590e-02] + # [ inf -1.13090565e-01 -4.93936523e-02 -2.74040590e-02 -1.57316125e-03]] + # h_s = + # [[ 0.00000000e+00 -8.62422934e+14 -3.76136070e+14 -1.61053099e+15 inf] + # [ -8.87390092e+14 2.59936788e-01 1.22409408e-01 7.97392657e-02 6.16999831e-02] + # [ -3.87924869e+14 1.46564082e-01 7.02793126e-02 4.62059958e-02 3.58607610e-02] + # [ -1.66844106e+15 1.18356989e-01 5.82970744e-02 3.92688958e-02 3.13685586e-02] + # [ inf 1.25606419e-01 6.35218804e-02 4.41902963e-02 3.69195895e-02]] + # k = find(w_1==w_2) # h_d(k) = h_dii @@ -915,7 +1059,6 @@ class SpecData1D(PlotData): # by pab 13.08.2002 - # TODO % Replace inputs with options structure # TODO % Can be improved further. method = 'apstochastic' @@ -3143,12 +3286,12 @@ class SpecData1D(PlotData): Parameters ---------- - nr : int + nr: int order of moments (recomended maximum 4) even : bool False for all moments, True for only even orders - j : int + j: int 0 or 1 Returns diff --git a/wafo/spectrum/models.py b/wafo/spectrum/models.py index 8799772..d753ecd 100644 --- a/wafo/spectrum/models.py +++ b/wafo/spectrum/models.py @@ -497,7 +497,7 @@ class Jonswap(ModelSpectrum): >>> import wafo.spectrum.models as wsm >>> S = wsm.Jonswap(Hm0=7, Tp=11,gamma=1) >>> S2 = wsm.Bretschneider(Hm0=7, Tp=11) - >>> w = plb.linspace(0,5) + >>> w = np.linspace(0,5) >>> all(np.abs(S(w)-S2(w))<1.e-7) True diff --git a/wafo/wave_theory/dispersion_relation.py b/wafo/wave_theory/dispersion_relation.py index b6eca61..818aa8c 100644 --- a/wafo/wave_theory/dispersion_relation.py +++ b/wafo/wave_theory/dispersion_relation.py @@ -9,7 +9,7 @@ import numpy as np from wafo.misc import lazywhere from numpy import (atleast_1d, sqrt, ones_like, zeros_like, arctan2, where, tanh, sin, cos, sign, inf, - flatnonzero, finfo, cosh) + flatnonzero, cosh) __all__ = ['k2w', 'w2k'] @@ -25,8 +25,7 @@ def _assert_warn(cond, msg): def k2w(k1, k2=0e0, h=inf, g=9.81, u1=0e0, u2=0e0): - """ Translates from wave number to frequency - using the dispersion relation + """ Translates from wave number to frequency using the dispersion relation Parameters ---------- @@ -73,9 +72,9 @@ def k2w(k1, k2=0e0, h=inf, g=9.81, u1=0e0, u2=0e0): array([ 0.13914927, 1.43498213, 2.00551724]) """ - k1i, k2i, hi, gi, u1i, u2i = atleast_1d(k1, k2, h, g, u1, u2) + k1i, k2i, hi, gi, u1i, u2i = np.broadcast_arrays(k1, k2, h, g, u1, u2) - if k1i.size == 0: + if np.size(k1i) == 0: return zeros_like(k1i) ku1 = k1i * u1i ku2 = k2i * u2i @@ -96,10 +95,9 @@ def k2w(k1, k2=0e0, h=inf, g=9.81, u1=0e0, u2=0e0): return w, theta -def w2k(w, theta=0.0, h=inf, g=9.81, count_limit=100): +def w2k(w, theta=0.0, h=inf, g=9.81, count_limit=100, rtol=1e-7, atol=1e-14): """ - Translates from frequency to wave number - using the dispersion relation + Translates from frequency to wave number using the dispersion relation Parameters ---------- @@ -145,22 +143,19 @@ def w2k(w, theta=0.0, h=inf, g=9.81, count_limit=100): -------- k2w """ - wi, th, hi, gi = atleast_1d(w, theta, h, g) - + gi = atleast_1d(g) + wi, th, hi = np.broadcast_arrays(w, theta, h) if wi.size == 0: return zeros_like(wi) k = 1.0 * sign(wi) * wi ** 2.0 / gi[0] # deep water - if (hi > 10. ** 25).all(): + if (hi > 1e25).all(): k2 = k * sin(th) * gi[0] / gi[-1] # size np x nf k1 = k * cos(th) return k1, k2 _assert(gi.size == 1, 'Finite depth in combination with 3D normalization' ' (len(g)=2) is not implemented yet.') - find = flatnonzero - eps = finfo(float).eps - oshape = k.shape wi, k, hi = wi.ravel(), k.ravel(), hi.ravel() @@ -168,11 +163,11 @@ def w2k(w, theta=0.0, h=inf, g=9.81, count_limit=100): # Permit no more than count_limit iterations. hi = hi * ones_like(k) hn = zeros_like(k) - ix = find((wi < 0) | (0 < wi)) + ix = flatnonzero(((wi < 0) | (0 < wi)) & (hi < 1e25)) # Break out of the iteration loop for three reasons: - # 1) the last update is very small (compared to x) - # 2) the last update is very small (compared to sqrt(eps)) + # 1) the last update is very small (compared to k*rtol) + # 2) the last update is very small (compared to atol) # 3) There are more than 100 iterations. This should NEVER happen. count = 0 while (ix.size > 0 and count < count_limit): @@ -180,13 +175,12 @@ def w2k(w, theta=0.0, h=inf, g=9.81, count_limit=100): kh = ki * hi[ix] coshkh2 = lazywhere(np.abs(kh) < 350, (kh, ), lambda kh: cosh(kh) ** 2.0, fillvalue=np.inf) - hn[ix] = (ki * tanh(kh) - wi[ix] ** 2.0 / gi) / \ - (tanh(kh) + kh / coshkh2) + hn[ix] = (ki * tanh(kh) - wi[ix] ** 2.0 / gi) / (tanh(kh) + kh / coshkh2) knew = ki - hn[ix] # Make sure that the current guess is not zero. # When Newton's Method suggests steps that lead to zero guesses # take a step 9/10ths of the way to zero: - ksmall = find(np.abs(knew) == 0) + ksmall = flatnonzero((np.abs(knew) == 0) | (np.isnan(knew)) ) if ksmall.size > 0: knew[ksmall] = ki[ksmall] / 10.0 hn[ix[ksmall]] = ki[ksmall] - knew[ksmall] @@ -195,8 +189,8 @@ def w2k(w, theta=0.0, h=inf, g=9.81, count_limit=100): # disp(['Iteration ',num2str(count),' Number of points left: ' # num2str(length(ix)) ]), - ix = find((np.abs(hn) > sqrt(eps) * np.abs(k)) * - np.abs(hn) > sqrt(eps)) + ix = flatnonzero((np.abs(hn) > rtol * np.abs(k)) * + (np.abs(hn) > atol)) count += 1 max_err = np.max(hn[ix]) if np.any(ix) else 0 _assert_warn(count < count_limit, 'W2K did not converge. ' diff --git a/wafo/wave_theory/tests/test_dispersion_relation.py b/wafo/wave_theory/tests/test_dispersion_relation.py index 4a92fa3..c76b4d5 100644 --- a/wafo/wave_theory/tests/test_dispersion_relation.py +++ b/wafo/wave_theory/tests/test_dispersion_relation.py @@ -4,31 +4,54 @@ Created on 19. juli 2010 @author: pab ''' import numpy as np +from numpy.testing import assert_allclose from wafo.wave_theory.dispersion_relation import w2k, k2w # @UnusedImport def test_k2w_infinite_water_depth(): vals = k2w(np.arange(0.01, .5, 0.2))[0] - true_vals = np.array([0.3132092, 1.43530485, 2.00551739]) - assert((np.abs(vals - true_vals) < 1e-7).all()) + true_vals = np.array([0.3132092, 1.43530485, 2.00551739]) + assert_allclose(vals, true_vals) def test_k2w_finite_water_depth(): - vals = k2w(np.arange(0.01, .5, 0.2), h=20)[0] - true_vals = np.array([0.13914927, 1.43498213, 2.00551724]) - assert((np.abs(vals - true_vals) < 1e-7).all()) + vals, theta = k2w(np.arange(0.01, .5, 0.2), h=20) + true_vals = (0.13914927, 1.43498213, 2.00551724) + assert_allclose(vals, true_vals) + assert_allclose(theta, 0) + + +def test_k2w_finite_water_depth_with_negative_k(): + vals, theta = k2w(-np.arange(0.01, .5, 0.2), h=20) + true_vals = [0.13914927, 1.43498213, 2.00551724] + assert_allclose(vals, true_vals) + assert_allclose(theta, np.pi) def test_w2k_infinite_water_depth(): - vals = w2k(range(4))[0] - true_vals = np.array([0., 0.1019368, 0.4077472, 0.91743119]) - assert((np.abs(vals - true_vals) < 1e-7).all()) + vals, k2 = w2k(range(4)) + true_vals = np.array([0., 0.1019368, 0.4077472, 0.91743119]) + assert_allclose(vals, true_vals) + assert_allclose(k2, 0) +def test_w2k_infinite_water_depth_with_negative_w(): + vals, k2 = w2k(-np.arange(4)) + true_vals = -1 * np.array([0., 0.1019368, 0.4077472, 0.91743119]) + assert_allclose(vals, true_vals) + assert_allclose(k2, 0) def test_w2k_finite_water_depth(): - vals = w2k(range(4), h=20)[0] - true_vals = np.array([0., 0.10503601, 0.40774726, 0.91743119]) - assert((np.abs(vals - true_vals) < 1e-7).all()) + vals, k2 = w2k(range(4), h=20) + true_vals = np.array([0., 0.10503601, 0.40774726, 0.91743119]) + assert_allclose(vals, true_vals) + assert_allclose(k2, 0) + +def test_w2k_finite_water_depth_with_negative_w(): + vals, k2 = w2k(-np.arange(4), h=20) + true_vals = -1 * np.array([0., 0.10503601, 0.40774726, 0.91743119]) + assert_allclose(vals, true_vals) + assert_allclose(k2, 0) + if __name__ == '__main__': import nose