source: flex_extract.git/Documentation/html/_modules/disaggregation.html @ 6931f61

ctbtodev
Last change on this file since 6931f61 was 6931f61, checked in by anphi <anne.philipp@…>, 4 years ago

Update Onlinedocumentation after review of language editing
  • Property mode set to 100644
File size: 67.9 KB
Line 
1
2
3<!DOCTYPE html>
4<!--[if IE 8]><html class="no-js lt-ie9" lang="en" > <![endif]-->
5<!--[if gt IE 8]><!--> <html class="no-js" lang="en" > <!--<![endif]-->
6<head>
7  <meta charset="utf-8">
8 
9  <meta name="viewport" content="width=device-width, initial-scale=1.0">
10 
11  <title>disaggregation &mdash; flex_extract 7.1.2 documentation</title>
12 
13
14 
15 
16 
17 
18
19 
20  <script type="text/javascript" src="../_static/js/modernizr.min.js"></script>
21 
22   
23      <script type="text/javascript" id="documentation_options" data-url_root="../" src="../_static/documentation_options.js"></script>
24        <script src="../_static/jquery.js"></script>
25        <script src="../_static/underscore.js"></script>
26        <script src="../_static/doctools.js"></script>
27        <script src="../_static/language_data.js"></script>
28        <script async="async" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/latest.js?config=TeX-AMS-MML_HTMLorMML"></script>
29   
30    <script type="text/javascript" src="../_static/js/theme.js"></script>
31
32   
33
34 
35  <link rel="stylesheet" href="../_static/css/theme.css" type="text/css" />
36  <link rel="stylesheet" href="../_static/pygments.css" type="text/css" />
37  <link rel="stylesheet" href="../_static/css/custom.css" type="text/css" />
38  <link rel="stylesheet" href="../_static/css/theme_overrides.css" type="text/css" />
39    <link rel="index" title="Index" href="../genindex.html" />
40    <link rel="search" title="Search" href="../search.html" /> 
41</head>
42
43<body class="wy-body-for-nav">
44
45   
46  <div class="wy-grid-for-nav">
47   
48    <nav data-toggle="wy-nav-shift" class="wy-nav-side">
49      <div class="wy-side-scroll">
50        <div class="wy-side-nav-search" >
51         
52
53         
54            <a href="../index.html" class="icon icon-home"> flex_extract
55         
56
57         
58          </a>
59
60         
61           
62           
63              <div class="version">
64                7.1.2
65              </div>
66           
67         
68
69         
70<div role="search">
71  <form id="rtd-search-form" class="wy-form" action="../search.html" method="get">
72    <input type="text" name="q" placeholder="Search docs" />
73    <input type="hidden" name="check_keywords" value="yes" />
74    <input type="hidden" name="area" value="default" />
75  </form>
76</div>
77
78         
79        </div>
80
81        <div class="wy-menu wy-menu-vertical" data-spy="affix" role="navigation" aria-label="main navigation">
82         
83           
84           
85             
86           
87           
88              <p class="caption"><span class="caption-text">Table of Contents:</span></p>
89<ul>
90<li class="toctree-l1"><a class="reference internal" href="../ecmwf_data.html">ECMWF Data</a></li>
91<li class="toctree-l1"><a class="reference internal" href="../installation.html">Installation</a></li>
92<li class="toctree-l1"><a class="reference internal" href="../quick_start.html">Usage</a></li>
93<li class="toctree-l1"><a class="reference internal" href="../documentation.html">Code-Level Documentation</a></li>
94<li class="toctree-l1"><a class="reference internal" href="../evaluation.html">Evaluation</a></li>
95<li class="toctree-l1"><a class="reference internal" href="../dev_guide.html">Developer Guide</a></li>
96<li class="toctree-l1"><a class="reference internal" href="../changelog.html">Changelog</a></li>
97<li class="toctree-l1"><a class="reference internal" href="../support.html">Support</a></li>
98<li class="toctree-l1"><a class="reference internal" href="../Support/faq.html">FAQ - Frequently asked questions</a></li>
99<li class="toctree-l1"><a class="reference internal" href="../authors.html">Developer Team</a></li>
100</ul>
101
102           
103         
104        </div>
105      </div>
106    </nav>
107
108    <section data-toggle="wy-nav-shift" class="wy-nav-content-wrap">
109
110     
111      <nav class="wy-nav-top" aria-label="top navigation">
112       
113          <i data-toggle="wy-nav-top" class="fa fa-bars"></i>
114          <a href="../index.html">flex_extract</a>
115       
116      </nav>
117
118
119      <div class="wy-nav-content">
120       
121        <div class="rst-content">
122       
123         
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139<div role="navigation" aria-label="breadcrumbs navigation">
140
141  <ul class="wy-breadcrumbs">
142   
143      <li><a href="../index.html">Docs</a> &raquo;</li>
144       
145          <li><a href="index.html">Module code</a> &raquo;</li>
146       
147      <li>disaggregation</li>
148   
149   
150      <li class="wy-breadcrumbs-aside">
151       
152      </li>
153   
154  </ul>
155
156 
157  <hr/>
158</div>
159          <div role="main" class="document" itemscope="itemscope" itemtype="http://schema.org/Article">
160           <div itemprop="articleBody">
161           
162  <h1>Source code for disaggregation</h1><div class="highlight"><pre>
163<span></span><span class="ch">#!/usr/bin/env python3</span>
164<span class="c1"># -*- coding: utf-8 -*-</span>
165<span class="c1">#*******************************************************************************</span>
166<span class="c1"># @Author: Anne Philipp (University of Vienna)</span>
167<span class="c1">#</span>
168<span class="c1"># @Date: March 2018</span>
169<span class="c1">#</span>
170<span class="c1"># @Change History:</span>
171<span class="c1">#</span>
172<span class="c1">#    November 2015 - Leopold Haimberger (University of Vienna):</span>
173<span class="c1">#        - migration of the methods dapoly and darain from Fortran</span>
174<span class="c1">#          (flex_extract_v6 and earlier) to Python</span>
175<span class="c1">#</span>
176<span class="c1">#    April 2018 - Anne Philipp (University of Vienna):</span>
177<span class="c1">#        - applied PEP8 style guide</span>
178<span class="c1">#        - added structured documentation</span>
179<span class="c1">#        - outsourced the disaggregation functions dapoly and darain</span>
180<span class="c1">#          to a new module named disaggregation</span>
181<span class="c1">#        - added the new disaggregation method for precipitation</span>
182<span class="c1">#</span>
183<span class="c1"># @License:</span>
184<span class="c1">#    (C) Copyright 2014-2019.</span>
185<span class="c1">#    Anne Philipp, Leopold Haimberger</span>
186<span class="c1">#</span>
187<span class="c1">#    SPDX-License-Identifier: CC-BY-4.0</span>
188<span class="c1">#</span>
189<span class="c1">#    This work is licensed under the Creative Commons Attribution 4.0</span>
190<span class="c1">#    International License. To view a copy of this license, visit</span>
191<span class="c1">#    http://creativecommons.org/licenses/by/4.0/ or send a letter to</span>
192<span class="c1">#    Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.</span>
193<span class="c1">#</span>
194<span class="c1"># @Methods:</span>
195<span class="c1">#    - dapoly</span>
196<span class="c1">#    - darain</span>
197<span class="c1">#    - IA3</span>
198<span class="c1">#*******************************************************************************</span>
199<span class="sd">&#39;&#39;&#39;Disaggregation of deaccumulated flux data from an ECMWF model FG field.</span>
200
201<span class="sd">Initially the flux data to be concerned are:</span>
202<span class="sd">    - large-scale precipitation</span>
203<span class="sd">    - convective precipitation</span>
204<span class="sd">    - surface sensible heat flux</span>
205<span class="sd">    - surface solar radiation</span>
206<span class="sd">    - u stress</span>
207<span class="sd">    - v stress</span>
208
209<span class="sd">Different versions of disaggregation is provided for rainfall</span>
210<span class="sd">data (darain, modified linear) and the surface fluxes and</span>
211<span class="sd">stress data (dapoly, cubic polynomial).</span>
212<span class="sd">&#39;&#39;&#39;</span>
213
214<span class="c1"># ------------------------------------------------------------------------------</span>
215<span class="c1"># MODULES</span>
216<span class="c1"># ------------------------------------------------------------------------------</span>
217
218<span class="c1"># ------------------------------------------------------------------------------</span>
219<span class="c1"># FUNCTIONS</span>
220<span class="c1"># ------------------------------------------------------------------------------</span>
221<div class="viewcode-block" id="dapoly"><a class="viewcode-back" href="../Documentation/Api/api_python.html#disaggregation.dapoly">[docs]</a><span class="k">def</span> <span class="nf">dapoly</span><span class="p">(</span><span class="n">alist</span><span class="p">):</span>
222    <span class="sd">&quot;&quot;&quot;Cubic polynomial interpolation of deaccumulated fluxes.</span>
223
224<span class="sd">    Interpolation of deaccumulated fluxes of an ECMWF model FG field</span>
225<span class="sd">    using a cubic polynomial solution which conserves the integrals</span>
226<span class="sd">    of the fluxes within each timespan.</span>
227<span class="sd">    Disaggregation is done for 4 accumluated timespans which</span>
228<span class="sd">    generates a new, disaggregated value which is output at the</span>
229<span class="sd">    central point of the 4 accumulation timespans.</span>
230<span class="sd">    This new point is used for linear interpolation of the complete</span>
231<span class="sd">    timeseries afterwards.</span>
232
233<span class="sd">    Parameters</span>
234<span class="sd">    ----------</span>
235<span class="sd">    alist : list of array of float</span>
236<span class="sd">        List of 4 timespans as 2-dimensional, horizontal fields.</span>
237<span class="sd">        E.g. [[array_t1], [array_t2], [array_t3], [array_t4]]</span>
238
239<span class="sd">    Return</span>
240<span class="sd">    ------</span>
241<span class="sd">    nfield : array of float</span>
242<span class="sd">        Interpolated flux at central point of accumulation timespan.</span>
243
244<span class="sd">    Note</span>
245<span class="sd">    ----</span>
246<span class="sd">    March 2000    : P. JAMES</span>
247<span class="sd">        Original author</span>
248
249<span class="sd">    June 2003     : A. BECK</span>
250<span class="sd">        Adaptations</span>
251
252<span class="sd">    November 2015 : Leopold Haimberger (University of Vienna)</span>
253<span class="sd">        Migration from Fortran to Python</span>
254
255<span class="sd">    &quot;&quot;&quot;</span>
256
257    <span class="n">pya</span> <span class="o">=</span> <span class="p">(</span><span class="n">alist</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">-</span> <span class="n">alist</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">+</span> <span class="mf">3.</span> <span class="o">*</span> <span class="p">(</span><span class="n">alist</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">alist</span><span class="p">[</span><span class="mi">2</span><span class="p">]))</span> <span class="o">/</span> <span class="mf">6.</span>
258    <span class="n">pyb</span> <span class="o">=</span> <span class="p">(</span><span class="n">alist</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span> <span class="o">+</span> <span class="n">alist</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span> <span class="o">/</span> <span class="mf">2.</span> <span class="o">-</span> <span class="n">alist</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="mf">9.</span> <span class="o">*</span> <span class="n">pya</span> <span class="o">/</span> <span class="mf">2.</span>
259    <span class="n">pyc</span> <span class="o">=</span> <span class="n">alist</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">alist</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">-</span> <span class="mf">7.</span> <span class="o">*</span> <span class="n">pya</span> <span class="o">/</span> <span class="mf">2.</span> <span class="o">-</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">pyb</span>
260    <span class="n">pyd</span> <span class="o">=</span> <span class="n">alist</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">-</span> <span class="n">pya</span> <span class="o">/</span> <span class="mf">4.</span> <span class="o">-</span> <span class="n">pyb</span> <span class="o">/</span> <span class="mf">3.</span> <span class="o">-</span> <span class="n">pyc</span> <span class="o">/</span> <span class="mf">2.</span>
261    <span class="n">nfield</span> <span class="o">=</span> <span class="mf">8.</span> <span class="o">*</span> <span class="n">pya</span> <span class="o">+</span> <span class="mf">4.</span> <span class="o">*</span> <span class="n">pyb</span> <span class="o">+</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">pyc</span> <span class="o">+</span> <span class="n">pyd</span>
262
263    <span class="k">return</span> <span class="n">nfield</span></div>
264
265
266<div class="viewcode-block" id="darain"><a class="viewcode-back" href="../Documentation/Api/api_python.html#disaggregation.darain">[docs]</a><span class="k">def</span> <span class="nf">darain</span><span class="p">(</span><span class="n">alist</span><span class="p">):</span>
267    <span class="sd">&quot;&quot;&quot;Linear interpolation of deaccumulated fluxes.</span>
268
269<span class="sd">    Interpolation of deaccumulated fluxes of an ECMWF model FG rainfall</span>
270<span class="sd">    field using a modified linear solution which conserves the integrals</span>
271<span class="sd">    of the fluxes within each timespan.</span>
272<span class="sd">    Disaggregation is done for 4 accumluated timespans which generates</span>
273<span class="sd">    a new, disaggregated value which is output at the central point</span>
274<span class="sd">    of the 4 accumulation timespans. This new point is used for linear</span>
275<span class="sd">    interpolation of the complete timeseries afterwards.</span>
276
277<span class="sd">    Parameters</span>
278<span class="sd">    ----------</span>
279<span class="sd">    alist : list of array of float</span>
280<span class="sd">        List of 4 timespans as 2-dimensional, horizontal fields.</span>
281<span class="sd">        E.g. [[array_t1], [array_t2], [array_t3], [array_t4]]</span>
282
283<span class="sd">    Return</span>
284<span class="sd">    ------</span>
285<span class="sd">    nfield : array of float</span>
286<span class="sd">        Interpolated flux at central point of accumulation timespan.</span>
287
288<span class="sd">    Note</span>
289<span class="sd">    ----</span>
290<span class="sd">    March 2000    : P. JAMES</span>
291<span class="sd">        Original author</span>
292
293<span class="sd">    June 2003     : A. BECK</span>
294<span class="sd">        Adaptations</span>
295
296<span class="sd">    November 2015 : Leopold Haimberger (University of Vienna)</span>
297<span class="sd">        Migration from Fortran to Python</span>
298<span class="sd">    &quot;&quot;&quot;</span>
299
300    <span class="n">xa</span> <span class="o">=</span> <span class="n">alist</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
301    <span class="n">xb</span> <span class="o">=</span> <span class="n">alist</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
302    <span class="n">xc</span> <span class="o">=</span> <span class="n">alist</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span>
303    <span class="n">xd</span> <span class="o">=</span> <span class="n">alist</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span>
304    <span class="n">xa</span><span class="p">[</span><span class="n">xa</span> <span class="o">&lt;</span> <span class="mf">0.</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.</span>
305    <span class="n">xb</span><span class="p">[</span><span class="n">xb</span> <span class="o">&lt;</span> <span class="mf">0.</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.</span>
306    <span class="n">xc</span><span class="p">[</span><span class="n">xc</span> <span class="o">&lt;</span> <span class="mf">0.</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.</span>
307    <span class="n">xd</span><span class="p">[</span><span class="n">xd</span> <span class="o">&lt;</span> <span class="mf">0.</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.</span>
308
309    <span class="n">xac</span> <span class="o">=</span> <span class="mf">0.5</span> <span class="o">*</span> <span class="n">xb</span>
310    <span class="n">mask</span> <span class="o">=</span> <span class="n">xa</span> <span class="o">+</span> <span class="n">xc</span> <span class="o">&gt;</span> <span class="mf">0.</span>
311    <span class="n">xac</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span> <span class="o">=</span> <span class="n">xb</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span> <span class="o">*</span> <span class="n">xc</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span> <span class="o">/</span> <span class="p">(</span><span class="n">xa</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span> <span class="o">+</span> <span class="n">xc</span><span class="p">[</span><span class="n">mask</span><span class="p">])</span>
312    <span class="n">xbd</span> <span class="o">=</span> <span class="mf">0.5</span> <span class="o">*</span> <span class="n">xc</span>
313    <span class="n">mask</span> <span class="o">=</span> <span class="n">xb</span> <span class="o">+</span> <span class="n">xd</span> <span class="o">&gt;</span> <span class="mf">0.</span>
314    <span class="n">xbd</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span> <span class="o">=</span> <span class="n">xb</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span> <span class="o">*</span> <span class="n">xc</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span> <span class="o">/</span> <span class="p">(</span><span class="n">xb</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span> <span class="o">+</span> <span class="n">xd</span><span class="p">[</span><span class="n">mask</span><span class="p">])</span>
315    <span class="n">nfield</span> <span class="o">=</span> <span class="n">xac</span> <span class="o">+</span> <span class="n">xbd</span>
316
317    <span class="k">return</span> <span class="n">nfield</span></div>
318
319<div class="viewcode-block" id="IA3"><a class="viewcode-back" href="../Documentation/Api/api_python.html#disaggregation.IA3">[docs]</a><span class="k">def</span> <span class="nf">IA3</span><span class="p">(</span><span class="n">g</span><span class="p">):</span>
320    <span class="sd">&quot;&quot;&quot; Interpolation with a non-negative geometric mean based algorithm.</span>
321
322<span class="sd">    The original grid is reconstructed by adding two sampling points in each</span>
323<span class="sd">    data series interval. This subgrid is used to keep all information during</span>
324<span class="sd">    the interpolation within the associated interval. Additionally, an advanced</span>
325<span class="sd">    monotonicity filter is applied to improve the monotonicity properties of</span>
326<span class="sd">    the series.</span>
327
328<span class="sd">    Note</span>
329<span class="sd">    ----</span>
330<span class="sd">    (C) Copyright 2017-2019</span>
331<span class="sd">    Sabine Hittmeir, Anne Philipp, Petra Seibert</span>
332
333<span class="sd">    This work is licensed under the Creative Commons Attribution 4.0</span>
334<span class="sd">    International License. To view a copy of this license, visit</span>
335<span class="sd">    http://creativecommons.org/licenses/by/4.0/ or send a letter to</span>
336<span class="sd">    Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.</span>
337
338<span class="sd">    Parameters</span>
339<span class="sd">    ----------</span>
340<span class="sd">    g : list of float</span>
341<span class="sd">        Complete data series that will be interpolated having</span>
342<span class="sd">        the dimension of the original raw series.</span>
343
344<span class="sd">    Return</span>
345<span class="sd">    ------</span>
346<span class="sd">    f : list of float</span>
347<span class="sd">        The interpolated data series with additional subgrid points.</span>
348<span class="sd">        Its dimension is equal to the length of the input data series</span>
349<span class="sd">        times three.</span>
350
351
352<span class="sd">    References</span>
353<span class="sd">    ----------</span>
354<span class="sd">    For more information see article:</span>
355<span class="sd">    Hittmeir, S.; Philipp, A.; Seibert, P. (2017): A conservative</span>
356<span class="sd">    interpolation scheme for extensive quantities with application to the</span>
357<span class="sd">    Lagrangian particle dispersion model FLEXPART.,</span>
358<span class="sd">    Geoscientific Model Development</span>
359<span class="sd">    &quot;&quot;&quot;</span>
360
361    <span class="c1">#######################  variable description #############################</span>
362    <span class="c1">#                                                                         #</span>
363    <span class="c1"># i      - index variable for looping over the data series                #</span>
364    <span class="c1"># g      - input data series                                              #</span>
365    <span class="c1"># f      - interpolated and filtered data series with additional          #</span>
366    <span class="c1">#          grid points                                                    #</span>
367    <span class="c1"># fi     - function value at position i, f_i                              #</span>
368    <span class="c1"># fi1    - first  sub-grid function value f_i^1                           #</span>
369    <span class="c1"># fi2    - second sub-grid function value f_i^2                           #</span>
370    <span class="c1"># fip1   - next function value at position i+1, f_(i+1)                   #</span>
371    <span class="c1"># dt     - time step                                                      #</span>
372    <span class="c1"># fmon   - monotonicity filter                                            #</span>
373    <span class="c1">#                                                                         #</span>
374    <span class="c1">###########################################################################</span>
375
376
377    <span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
378
379    <span class="c1"># time step</span>
380    <span class="n">dt</span> <span class="o">=</span> <span class="mf">1.0</span>
381
382    <span class="c1">############### Non-negative Geometric Mean Based Algorithm ###############</span>
383
384    <span class="c1"># for the left boundary the following boundary condition is valid:</span>
385    <span class="c1"># the value at t=0 of the interpolation algorithm coincides with the</span>
386    <span class="c1"># first data value according to the persistence hypothesis</span>
387    <span class="n">f</span> <span class="o">=</span> <span class="p">[</span><span class="n">g</span><span class="p">[</span><span class="mi">0</span><span class="p">]]</span>
388
389    <span class="c1"># compute two first sub-grid intervals without monotonicity check</span>
390    <span class="c1"># go through the data series and extend each interval by two sub-grid</span>
391    <span class="c1"># points and interpolate the corresponding data values</span>
392    <span class="c1"># except for the last interval due to boundary conditions</span>
393    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mi">2</span><span class="p">):</span>
394
395        <span class="c1"># as a requirement:</span>
396        <span class="c1"># if there is a zero data value such that g[i]=0, then the whole</span>
397        <span class="c1"># interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0</span>
398        <span class="c1"># according to Eq. (6)</span>
399        <span class="k">if</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">==</span> <span class="mf">0.</span><span class="p">:</span>
400            <span class="n">f</span><span class="o">.</span><span class="n">extend</span><span class="p">([</span><span class="mf">0.</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="mf">0.</span><span class="p">])</span>
401
402        <span class="c1"># otherwise the sub-grid values are calculated and added to the list</span>
403        <span class="k">else</span><span class="p">:</span>
404            <span class="c1"># temporal save of last value in interpolated list</span>
405            <span class="c1"># since it is the left boundary and hence the new (fi) value</span>
406            <span class="n">fi</span> <span class="o">=</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
407
408            <span class="c1"># the value at the end of the interval (fip1) is prescribed by the</span>
409            <span class="c1"># geometric mean, restricted such that non-negativity is guaranteed</span>
410            <span class="c1"># according to Eq. (25)</span>
411            <span class="n">fip1</span> <span class="o">=</span> <span class="nb">min</span><span class="p">(</span><span class="mf">3.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span><span class="p">],</span> <span class="mf">3.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span> <span class="o">+</span> <span class="mi">1</span><span class="p">],</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">g</span><span class="p">[</span><span class="n">i</span> <span class="o">+</span> <span class="mi">1</span><span class="p">]</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span><span class="p">]))</span>
412
413            <span class="c1"># the function value at the first sub-grid point (fi1) is determined</span>
414            <span class="c1"># according to the equal area condition with Eq. (19)</span>
415            <span class="n">fi1</span> <span class="o">=</span> <span class="mf">3.</span><span class="o">/</span><span class="mf">2.</span><span class="o">*</span><span class="n">g</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">-</span><span class="mf">5.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">fip1</span><span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">fi</span>
416
417            <span class="c1"># the function value at the second sub-grid point (fi2) is determined</span>
418            <span class="c1"># according Eq. (18)</span>
419            <span class="n">fi2</span> <span class="o">=</span> <span class="n">fi1</span><span class="o">+</span><span class="mf">1.</span><span class="o">/</span><span class="mf">3.</span><span class="o">*</span><span class="p">(</span><span class="n">fip1</span><span class="o">-</span><span class="n">fi</span><span class="p">)</span>
420
421            <span class="c1"># add next interval of interpolated (sub-)grid values</span>
422            <span class="n">f</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">fi1</span><span class="p">)</span>
423            <span class="n">f</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">fi2</span><span class="p">)</span>
424            <span class="n">f</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">fip1</span><span class="p">)</span>
425
426    <span class="c1"># compute rest of the data series intervals</span>
427    <span class="c1"># go through the data series and extend each interval by two sub-grid</span>
428    <span class="c1"># points and interpolate the corresponding data values</span>
429    <span class="c1"># except for the last interval due to boundary conditions</span>
430    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">2</span><span class="p">,</span> <span class="nb">len</span><span class="p">(</span><span class="n">g</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
431
432        <span class="c1"># as a requirement:</span>
433        <span class="c1"># if there is a zero data value such that g[i]=0, then the whole</span>
434        <span class="c1"># interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0</span>
435        <span class="c1"># according to Eq. (6)</span>
436        <span class="k">if</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">==</span> <span class="mf">0.</span><span class="p">:</span>
437            <span class="c1"># apply monotonicity filter for interval before</span>
438            <span class="c1"># check if there is &quot;M&quot; or &quot;W&quot; shape</span>
439            <span class="k">if</span>     <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">6</span><span class="p">])</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">])</span> <span class="o">==</span> <span class="o">-</span><span class="mi">1</span> \
440               <span class="ow">and</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">])</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">])</span> <span class="o">==</span> <span class="o">-</span><span class="mi">1</span> \
441               <span class="ow">and</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">])</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">])</span> <span class="o">==</span> <span class="o">-</span><span class="mi">1</span><span class="p">:</span>
442
443                <span class="c1"># the monotonicity filter corrects the value at (fim1) by</span>
444                <span class="c1"># substituting (fim1) with (fmon), see Eq. (27), (28) and (29)</span>
445                <span class="n">fmon</span> <span class="o">=</span> <span class="nb">min</span><span class="p">(</span><span class="mf">3.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span> <span class="o">-</span> <span class="mi">2</span><span class="p">],</span>
446                           <span class="mf">3.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span> <span class="o">-</span> <span class="mi">1</span><span class="p">],</span>
447                           <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="nb">max</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="p">(</span><span class="mf">18.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span> <span class="o">-</span> <span class="mi">2</span><span class="p">]</span> <span class="o">-</span> <span class="mf">5.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">7</span><span class="p">])</span> <span class="o">*</span>
448                                       <span class="p">(</span><span class="mf">18.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span> <span class="o">-</span> <span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="mf">5.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]))))</span>
449
450                <span class="c1"># recomputation of the sub-grid interval values while the</span>
451                <span class="c1"># interval boundaries (fi) and (fip2) remains unchanged</span>
452                <span class="c1"># see Eq. (18) and (19)</span>
453                <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">]</span> <span class="o">=</span> <span class="n">fmon</span>
454                <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">6</span><span class="p">]</span> <span class="o">=</span> <span class="mf">3.</span><span class="o">/</span><span class="mf">2.</span><span class="o">*</span><span class="n">g</span><span class="p">[</span><span class="n">i</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span><span class="o">-</span><span class="mf">5.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">fmon</span><span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">7</span><span class="p">]</span>
455                <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">]</span> <span class="o">=</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">6</span><span class="p">]</span><span class="o">+</span><span class="p">(</span><span class="n">fmon</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">7</span><span class="p">])</span><span class="o">/</span><span class="mf">3.</span>
456                <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="mf">3.</span><span class="o">/</span><span class="mf">2.</span><span class="o">*</span><span class="n">g</span><span class="p">[</span><span class="n">i</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="mf">5.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">fmon</span>
457                <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span><span class="o">+</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">fmon</span><span class="p">)</span><span class="o">/</span><span class="mf">3.</span>
458
459            <span class="n">f</span><span class="o">.</span><span class="n">extend</span><span class="p">([</span><span class="mf">0.</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="mf">0.</span><span class="p">])</span>
460
461        <span class="c1"># otherwise the sub-grid values are calculated and added to the list</span>
462        <span class="k">else</span><span class="p">:</span>
463            <span class="c1"># temporal save of last value in interpolated list</span>
464            <span class="c1"># since it is the left boundary and hence the new (fi) value</span>
465            <span class="n">fi</span> <span class="o">=</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
466
467            <span class="c1"># the value at the end of the interval (fip1) is prescribed by the</span>
468            <span class="c1"># geometric mean, restricted such that non-negativity is guaranteed</span>
469            <span class="c1"># according to Eq. (25)</span>
470            <span class="n">fip1</span> <span class="o">=</span> <span class="nb">min</span><span class="p">(</span><span class="mf">3.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span><span class="p">],</span> <span class="mf">3.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span> <span class="o">+</span> <span class="mi">1</span><span class="p">],</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">g</span><span class="p">[</span><span class="n">i</span> <span class="o">+</span> <span class="mi">1</span><span class="p">]</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span><span class="p">]))</span>
471
472            <span class="c1"># the function value at the first sub-grid point (fi1) is determined</span>
473            <span class="c1"># according to the equal area condition with Eq. (19)</span>
474            <span class="n">fi1</span> <span class="o">=</span> <span class="mf">3.</span><span class="o">/</span><span class="mf">2.</span><span class="o">*</span><span class="n">g</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">-</span><span class="mf">5.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">fip1</span><span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">fi</span>
475
476            <span class="c1"># the function value at the second sub-grid point (fi2) is determined</span>
477            <span class="c1"># according Eq. (18)</span>
478            <span class="n">fi2</span> <span class="o">=</span> <span class="n">fi1</span><span class="o">+</span><span class="mf">1.</span><span class="o">/</span><span class="mf">3.</span><span class="o">*</span><span class="p">(</span><span class="n">fip1</span><span class="o">-</span><span class="n">fi</span><span class="p">)</span>
479
480            <span class="c1"># apply monotonicity filter for interval before</span>
481            <span class="c1"># check if there is &quot;M&quot; or &quot;W&quot; shape</span>
482            <span class="k">if</span>     <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">6</span><span class="p">])</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">])</span> <span class="o">==</span> <span class="o">-</span><span class="mi">1</span> \
483               <span class="ow">and</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">])</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">])</span> <span class="o">==</span> <span class="o">-</span><span class="mi">1</span> \
484               <span class="ow">and</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">])</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">])</span> <span class="o">==</span> <span class="o">-</span><span class="mi">1</span><span class="p">:</span>
485
486                <span class="c1"># the monotonicity filter corrects the value at (fim1) by</span>
487                <span class="c1"># substituting (fim1) with fmon, see Eq. (27), (28) and (29)</span>
488                <span class="n">fmon</span> <span class="o">=</span> <span class="nb">min</span><span class="p">(</span><span class="mf">3.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span> <span class="o">-</span> <span class="mi">2</span><span class="p">],</span>
489                           <span class="mf">3.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span> <span class="o">-</span> <span class="mi">1</span><span class="p">],</span>
490                           <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="nb">max</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="p">(</span><span class="mf">18.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span> <span class="o">-</span> <span class="mi">2</span><span class="p">]</span> <span class="o">-</span> <span class="mf">5.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">7</span><span class="p">])</span> <span class="o">*</span>
491                                       <span class="p">(</span><span class="mf">18.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="n">i</span> <span class="o">-</span> <span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="mf">5.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]))))</span>
492
493                <span class="c1"># recomputation of the sub-grid interval values while the</span>
494                <span class="c1"># interval boundaries (fi) and (fip2) remains unchanged</span>
495                <span class="c1"># see Eq. (18) and (19)</span>
496                <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">]</span> <span class="o">=</span> <span class="n">fmon</span>
497                <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">6</span><span class="p">]</span> <span class="o">=</span> <span class="mf">3.</span><span class="o">/</span><span class="mf">2.</span><span class="o">*</span><span class="n">g</span><span class="p">[</span><span class="n">i</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span><span class="o">-</span><span class="mf">5.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">fmon</span><span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">7</span><span class="p">]</span>
498                <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">]</span> <span class="o">=</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">6</span><span class="p">]</span><span class="o">+</span><span class="p">(</span><span class="n">fmon</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">7</span><span class="p">])</span><span class="o">/</span><span class="mf">3.</span>
499                <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="mf">3.</span><span class="o">/</span><span class="mf">2.</span><span class="o">*</span><span class="n">g</span><span class="p">[</span><span class="n">i</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="mf">5.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">fmon</span>
500                <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span><span class="o">+</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">fmon</span><span class="p">)</span><span class="o">/</span><span class="mf">3.</span>
501
502            <span class="c1"># add next interval of interpolated (sub-)grid values</span>
503            <span class="n">f</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">fi1</span><span class="p">)</span>
504            <span class="n">f</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">fi2</span><span class="p">)</span>
505            <span class="n">f</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">fip1</span><span class="p">)</span>
506
507    <span class="c1"># separate treatment of the final interval</span>
508
509    <span class="c1"># as a requirement:</span>
510    <span class="c1"># if there is a zero data value such that g[i]=0, then the whole</span>
511    <span class="c1"># interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0</span>
512    <span class="c1"># according to Eq. (6)</span>
513    <span class="k">if</span> <span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">==</span> <span class="mf">0.</span><span class="p">:</span>
514        <span class="c1"># apply monotonicity filter for interval before</span>
515        <span class="c1"># check if there is &quot;M&quot; or &quot;W&quot; shape</span>
516        <span class="k">if</span>     <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">6</span><span class="p">])</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">])</span> <span class="o">==</span> <span class="o">-</span><span class="mi">1</span> \
517           <span class="ow">and</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">])</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">])</span> <span class="o">==</span> <span class="o">-</span><span class="mi">1</span> \
518           <span class="ow">and</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">])</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">])</span> <span class="o">==</span> <span class="o">-</span><span class="mi">1</span><span class="p">:</span>
519
520            <span class="c1"># the monotonicity filter corrects the value at (fim1) by</span>
521            <span class="c1"># substituting (fim1) with (fmon), see Eq. (27), (28) and (29)</span>
522            <span class="n">fmon</span> <span class="o">=</span> <span class="nb">min</span><span class="p">(</span><span class="mf">3.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">],</span>
523                       <span class="mf">3.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">],</span>
524                       <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="nb">max</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="p">(</span><span class="mf">18.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span> <span class="o">-</span> <span class="mf">5.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">7</span><span class="p">])</span> <span class="o">*</span>
525                                   <span class="p">(</span><span class="mf">18.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span> <span class="o">-</span> <span class="mf">5.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]))))</span>
526
527            <span class="c1"># recomputation of the sub-grid interval values while the</span>
528            <span class="c1"># interval boundaries (fi) and (fip2) remains unchanged</span>
529            <span class="c1"># see Eq. (18) and (19)</span>
530            <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">]</span> <span class="o">=</span> <span class="n">fmon</span>
531            <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">6</span><span class="p">]</span> <span class="o">=</span> <span class="mf">3.</span><span class="o">/</span><span class="mf">2.</span><span class="o">*</span><span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span><span class="o">-</span><span class="mf">5.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">fmon</span><span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">7</span><span class="p">]</span>
532            <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">]</span> <span class="o">=</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">6</span><span class="p">]</span><span class="o">+</span><span class="p">(</span><span class="n">fmon</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">7</span><span class="p">])</span><span class="o">/</span><span class="mf">3.</span>
533            <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="mf">3.</span><span class="o">/</span><span class="mf">2.</span><span class="o">*</span><span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span><span class="o">-</span><span class="mf">5.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">fmon</span>
534            <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span><span class="o">+</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">fmon</span><span class="p">)</span><span class="o">/</span><span class="mf">3.</span>
535
536        <span class="n">f</span><span class="o">.</span><span class="n">extend</span><span class="p">([</span><span class="mf">0.</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="mf">0.</span><span class="p">])</span>
537
538    <span class="c1"># otherwise the sub-grid values are calculated and added to the list</span>
539    <span class="c1"># using the persistence hypothesis as boundary condition</span>
540    <span class="k">else</span><span class="p">:</span>
541        <span class="c1"># temporal save of last value in interpolated list</span>
542        <span class="c1"># since it is the left boundary and hence the new (fi) value</span>
543        <span class="n">fi</span> <span class="o">=</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
544        <span class="c1"># since last interval in series, last value is also fip1</span>
545        <span class="n">fip1</span> <span class="o">=</span> <span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
546        <span class="c1"># the function value at the first sub-grid point (fi1) is determined</span>
547        <span class="c1"># according to the equal area condition with Eq. (19)</span>
548        <span class="n">fi1</span> <span class="o">=</span> <span class="mf">3.</span><span class="o">/</span><span class="mf">2.</span><span class="o">*</span><span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="mf">5.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">fip1</span><span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">fi</span>
549        <span class="c1"># the function value at the second sub-grid point (fi2) is determined</span>
550        <span class="c1"># according Eq. (18)</span>
551        <span class="n">fi2</span> <span class="o">=</span> <span class="n">fi1</span><span class="o">+</span><span class="n">dt</span><span class="o">/</span><span class="mf">3.</span><span class="o">*</span><span class="p">(</span><span class="n">fip1</span><span class="o">-</span><span class="n">fi</span><span class="p">)</span>
552
553        <span class="c1"># apply monotonicity filter for interval before</span>
554        <span class="c1"># check if there is &quot;M&quot; or &quot;W&quot; shape</span>
555        <span class="k">if</span>     <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">6</span><span class="p">])</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">])</span> <span class="o">==</span> <span class="o">-</span><span class="mi">1</span> \
556           <span class="ow">and</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">])</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">])</span> <span class="o">==</span> <span class="o">-</span><span class="mi">1</span> \
557           <span class="ow">and</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">])</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">])</span> <span class="o">==</span> <span class="o">-</span><span class="mi">1</span><span class="p">:</span>
558
559            <span class="c1"># the monotonicity filter corrects the value at (fim1) by</span>
560            <span class="c1"># substituting (fim1) with (fmon), see Eq. (27), (28) and (29)</span>
561            <span class="n">fmon</span> <span class="o">=</span> <span class="nb">min</span><span class="p">(</span><span class="mf">3.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">],</span>
562                       <span class="mf">3.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">],</span>
563                       <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="nb">max</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="p">(</span><span class="mf">18.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span> <span class="o">-</span> <span class="mf">5.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">7</span><span class="p">])</span> <span class="o">*</span>
564                                   <span class="p">(</span><span class="mf">18.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span> <span class="o">-</span> <span class="mf">5.</span> <span class="o">/</span> <span class="mf">13.</span> <span class="o">*</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]))))</span>
565
566            <span class="c1"># recomputation of the sub-grid interval values while the</span>
567            <span class="c1"># interval boundaries (fi) and (fip2) remains unchanged</span>
568            <span class="c1"># see Eq. (18) and (19)</span>
569            <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">4</span><span class="p">]</span> <span class="o">=</span> <span class="n">fmon</span>
570            <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">6</span><span class="p">]</span> <span class="o">=</span> <span class="mf">3.</span><span class="o">/</span><span class="mf">2.</span><span class="o">*</span><span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span><span class="o">-</span><span class="mf">5.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">fmon</span><span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">7</span><span class="p">]</span>
571            <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">5</span><span class="p">]</span> <span class="o">=</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">6</span><span class="p">]</span><span class="o">+</span><span class="p">(</span><span class="n">fmon</span><span class="o">-</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">7</span><span class="p">])</span><span class="o">/</span><span class="mf">3.</span>
572            <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="mf">3.</span><span class="o">/</span><span class="mf">2.</span><span class="o">*</span><span class="n">g</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span><span class="o">-</span><span class="mf">5.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="mf">12.</span><span class="o">*</span><span class="n">fmon</span>
573            <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">3</span><span class="p">]</span><span class="o">+</span><span class="p">(</span><span class="n">f</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">fmon</span><span class="p">)</span><span class="o">/</span><span class="mf">3.</span>
574
575        <span class="c1"># add next interval of interpolated (sub-)grid values</span>
576        <span class="n">f</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">fi1</span><span class="p">)</span>
577        <span class="n">f</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">fi2</span><span class="p">)</span>
578        <span class="n">f</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">fip1</span><span class="p">)</span>
579
580    <span class="k">return</span> <span class="n">f</span></div>
581</pre></div>
582
583           </div>
584           
585          </div>
586          <footer>
587 
588
589  <hr/>
590
591  <div role="contentinfo">
592    <p>
593        &copy; Copyright 2020, Anne Philipp, Leopold Haimberger and Petra Seibert
594
595    </p>
596  </div>
597  Built with <a href="http://sphinx-doc.org/">Sphinx</a> using a <a href="https://github.com/rtfd/sphinx_rtd_theme">theme</a> provided by <a href="https://readthedocs.org">Read the Docs</a>.
598
599</footer>
600
601        </div>
602      </div>
603
604    </section>
605
606  </div>
607 
608
609
610  <script type="text/javascript">
611      jQuery(function () {
612          SphinxRtdTheme.Navigation.enable(true);
613      });
614  </script>
615
616 
617 
618   
619   
620
621</body>
622</html>
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG