diff --git a/.github/workflows/continuous-integration-workflow.yaml b/.github/workflows/continuous-integration-workflow.yaml new file mode 100644 index 0000000..f0937a2 --- /dev/null +++ b/.github/workflows/continuous-integration-workflow.yaml @@ -0,0 +1,37 @@ +name: Python System Tests with Coverage + +on: + push: + branches: + - main + - develop + pull_request: + branches: + - main + - develop + +jobs: + test: + runs-on: ubuntu-latest + + steps: + - name: Checkout code + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: '3.11, 3.12, 3.13' + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install .[dev] + + - name: Install FLOWERS + run: | + pip install .[dev] + + - name: Run system tests with coverage + run: | + pytest --cov=flowers --cov-fail-under=50 test/system \ No newline at end of file diff --git a/.gitignore b/.gitignore index 53ca5d6..bac859b 100644 --- a/.gitignore +++ b/.gitignore @@ -1,27 +1,2 @@ /.DS_Store /__pycache__ -/testing.py -/floris.mp4 -/flowers.mp4 -/project_solutions -/output -/solutions -floris_1.mp4 -floris_9.mp4 -flowers_1.mp4 -flowers_9.mp4 -array.sh -array_floris.sh -array_flowers.sh -gauss_floris.sh -gauss_flowers.sh -scaling_floris.sh -scaling_flowers.sh -debug.sh -case48_initial.csv -case48_flowers.csv -case48_floris.csv -/old_files -/figures -schematics.py -/.vscode \ No newline at end of file diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..4c8c37f --- /dev/null +++ b/LICENSE @@ -0,0 +1,199 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..137994a --- /dev/null +++ b/README.md @@ -0,0 +1,3 @@ +# FLOWERS - FLOW Estimation and Rose Superposition + +tbd \ No newline at end of file diff --git a/aep_analysis.py b/aep_analysis.py deleted file mode 100644 index c885f6a..0000000 --- a/aep_analysis.py +++ /dev/null @@ -1,306 +0,0 @@ -import numpy as np -import pickle -from shapely.geometry import Polygon, Point, LineString -import warnings - -import model_interface as inter -import tools as tl - -warnings.filterwarnings("ignore") - -""" -This file conducts three different numerical experiments to compare -AEP estimates and their sensitivity between the FLOWERS and Conventional -AEP models - -""" - -case = 0 - - -########################################################################### -# RANDOMIZED LAYOUT / WIND ROSE -########################################################################### -if case == 0: - # Analysis options - N_max = 200 - num_terms= 10 - wd_resolution = 5.0 - ws_avg = True - - # Initialization - file = 'solutions/aep' + str(case) + '.p' - wr_list = [] - N_turb = np.random.randint(2,501,N_max) - N_wr = np.random.randint(1,10,N_max) - - aep_flowers = np.zeros(N_max) - aep_floris = np.zeros(N_max) - time_flowers = np.zeros(N_max) - time_floris = np.zeros(N_max) - - # Gather wind roses - for i in np.arange(1,10): - wr_list.append(tl.load_wind_rose(i)) - - model = inter.AEPInterface(wr_list[0], np.zeros(2), np.zeros(2)) - - # Calculate AEP for random layout and wind rose - for idx in range(N_max): - print('Case ' + str(idx+1) + ' of ' + str(N_max)) - wind_rose = wr_list[N_wr[idx]-1] - layout_x, layout_y = tl.discrete_layout(n_turb=N_turb[idx], min_dist=3.0) - - model.reinitialize(wind_rose=wind_rose, layout_x=layout_x, layout_y=layout_y, num_terms=num_terms, wd_resolution=wd_resolution, ws_avg=ws_avg) - aep, time = model.compare_aep(display=False) - - aep_flowers[idx] = aep[0] - aep_floris[idx] = aep[1] - - time_flowers[idx] = time[0] - time_floris[idx] = time[1] - - var = (N_wr, N_turb) - - pickle.dump((var, aep_flowers, aep_floris, time_flowers, time_floris), open(file,'wb')) - -########################################################################### -# RANDOMIZED RESOLUTION -########################################################################### -if case == 1: - # Analysis options - N_max = 200 - N_turb = 50 - ws_avg = True - - # Initialization - file = 'solutions/aep' + str(case) + '.p' - wr_list = [] - N_wr = np.random.randint(1,10,N_max) - N_terms = np.random.randint(2,180,N_max) - N_bins = np.random.randint(4,360,N_max) - - aep_flowers = np.zeros(N_max) - aep_floris = np.zeros(N_max) - time_flowers = np.zeros(N_max) - time_floris = np.zeros(N_max) - - # Gather wind roses - for i in np.arange(1,10): - wr_list.append(tl.load_wind_rose(i)) - - model = inter.AEPInterface(wr_list[0], np.zeros(2), np.zeros(2)) - - # Calculate AEP for random layout, wind rose, and resolution - for idx in range(N_max): - print('Case ' + str(idx+1) + ' of ' + str(N_max)) - wind_rose = wr_list[N_wr[idx]-1] - layout_x, layout_y = tl.discrete_layout(n_turb=N_turb, min_dist=3.0) - - model.reinitialize(wind_rose=wind_rose, layout_x=layout_x, layout_y=layout_y, num_terms=N_terms[idx], wd_resolution=360/N_bins[idx], ws_avg=ws_avg) - aep, time = model.compare_aep(display=False) - - aep_flowers[idx] = aep[0] - aep_floris[idx] = aep[1] - - time_flowers[idx] = time[0] - time_floris[idx] = time[1] - - var = (N_terms, N_bins) - - pickle.dump((var, aep_flowers, aep_floris, time_flowers, time_floris), open(file,'wb')) - -########################################################################### -# LAYOUT MUTATION -########################################################################### -if case == 2: - # Analysis options - N = 25 - step = 126./5 - np.random.seed(4) - flowers_terms = [100,20,10] - conv_resolution = [1,5,10] - ws_avg = True - - # Define parameters - X = 126. * np.array([0.,0.,0.,7.,7.,7.,14.,14.,14.]) - Y = 126. * np.array([0.,7.,14.,0.,7.,14.,0.,7.,14.]) - wr = tl.load_wind_rose(7) - nt = len(X) - file = 'solutions/aep' + str(case) + '.p' - resolution = [flowers_terms, conv_resolution] - - # Store AEP - aep_flowers = np.ones((3,N+1)) - aep_park = np.ones((3,N+1)) - aep_gauss = np.ones((3,N+1)) - x_all = np.zeros((N+1,nt)) - y_all = np.zeros((N+1,nt)) - - # Define FLOWERS and Park models - model_high = inter.AEPInterface(wr, X, Y, conventional_model='park') - model_high.reinitialize(num_terms=flowers_terms[0], wd_resolution=conv_resolution[0], ws_avg=ws_avg) - - model_med = inter.AEPInterface(wr, X, Y, conventional_model='park') - model_med.reinitialize(num_terms=flowers_terms[1], wd_resolution=conv_resolution[1], ws_avg=ws_avg) - - model_low = inter.AEPInterface(wr, X, Y, conventional_model='park') - model_low.reinitialize(num_terms=flowers_terms[2], wd_resolution=conv_resolution[2], ws_avg=ws_avg) - - # Define Gauss models - gauss_high = inter.AEPInterface(wr, X, Y, conventional_model='gauss') - gauss_high.reinitialize(num_terms=flowers_terms[0], wd_resolution=conv_resolution[0], ws_avg=ws_avg) - - gauss_med = inter.AEPInterface(wr, X, Y, conventional_model='gauss') - gauss_med.reinitialize(num_terms=flowers_terms[1], wd_resolution=conv_resolution[1], ws_avg=ws_avg) - - gauss_low = inter.AEPInterface(wr, X, Y, conventional_model='gauss') - gauss_low.reinitialize(num_terms=flowers_terms[2], wd_resolution=conv_resolution[2], ws_avg=ws_avg) - - # Store initial data - aep_flowers[0,0] = model_high.compute_flowers_aep() - aep_flowers[1,0] = model_med.compute_flowers_aep() - aep_flowers[2,0] = model_low.compute_flowers_aep() - - aep_park[0,0] = model_high.compute_floris_aep() - aep_park[1,0] = model_med.compute_floris_aep() - aep_park[2,0] = model_low.compute_floris_aep() - - aep_gauss[0,0] = gauss_high.compute_floris_aep() - aep_gauss[1,0] = gauss_med.compute_floris_aep() - aep_gauss[2,0] = gauss_low.compute_floris_aep() - - x_all[0] = X - y_all[0] = Y - - for i in np.arange(1,N+1): - print('Case ' + str(i) + ' of ' + str(N)) - - # Add random Gaussian noise to turbine positions - X += np.random.normal(0.,step,nt) - Y += np.random.normal(0.,step,nt) - - model_high.reinitialize(layout_x=X, layout_y=Y) - model_med.reinitialize(layout_x=X, layout_y=Y) - model_low.reinitialize(layout_x=X, layout_y=Y) - gauss_high.reinitialize(layout_x=X, layout_y=Y) - gauss_med.reinitialize(layout_x=X, layout_y=Y) - gauss_low.reinitialize(layout_x=X, layout_y=Y) - - # Compute new AEP for each model - aep_flowers[0,i] = model_high.compute_flowers_aep() - aep_flowers[1,i] = model_med.compute_flowers_aep() - aep_flowers[2,i] = model_low.compute_flowers_aep() - - aep_park[0,i] = model_high.compute_floris_aep() - aep_park[1,i] = model_med.compute_floris_aep() - aep_park[2,i] = model_low.compute_floris_aep() - - aep_gauss[0,i] = gauss_high.compute_floris_aep() - aep_gauss[1,i] = gauss_med.compute_floris_aep() - aep_gauss[2,i] = gauss_low.compute_floris_aep() - - x_all[i] = np.copy(X) - y_all[i] = np.copy(Y) - - # Normalize AEP by initial value - aep_flowers /= np.expand_dims(aep_flowers[:,0],-1) - aep_park /= np.expand_dims(aep_park[:,0],-1) - aep_gauss /= np.expand_dims(aep_gauss[:,0],-1) - - pickle.dump((x_all, y_all, aep_flowers, aep_park, aep_gauss, resolution, wr), open(file,'wb')) - -########################################################################### -# SOLUTION SPACE SMOOTHNESS -########################################################################### -if case == 3: - # Resolution options - flowers_terms = [180,90,60,45,36,30,25,20,15,10,5] - conv_resolution = [1,2,3,4,5,6,8,10,12,15,18] - - # Load layout, boundaries, and rose - X0 = 126. * np.array([0.,0.,0.,7.,7.,7.,14.,14.,14.]) - Y0 = 126. * np.array([0.,7.,14.,0.,7.,14.,0.,7.,14.]) - idx = 4 - boundaries = [(0., 0.),(0, 14*126.),(14*126, 14*126.),(14*126, 0.)] - nx = 16*4+1 - ny = 16*4+1 - - wind_rose = tl.load_wind_rose(7) - - # Initialization - file = 'solutions/aep' + str(case) + '.p' - xmin = np.min(X0) - 126. - xmax = np.max(X0) + 126. - ymin = np.min(Y0) - 126. - ymax = np.max(Y0) + 126. - - model = inter.AEPInterface(wind_rose, X0, Y0, conventional_model='park') - gauss = inter.AEPInterface(wind_rose, X0, Y0, num_terms=3, conventional_model='gauss') - - xx = np.linspace(xmin,xmax,nx,endpoint=True) - yy = np.linspace(ymin,ymax,ny,endpoint=True) - xx,yy = np.meshgrid(xx,yy) - - aep_flowers = np.zeros((len(flowers_terms),nx,ny)) - aep_park = np.zeros((len(conv_resolution),nx,ny)) - aep_gauss = np.zeros((len(conv_resolution),nx,ny)) - - poly = Polygon(boundaries) - line = LineString(boundaries) - - for n in range(len(flowers_terms)): - print('Resolution ' + str(n+1) + ' of ' + str(len(flowers_terms))) - infeasible = np.zeros_like(xx) - - # Set resolution for models - model.reinitialize(num_terms=flowers_terms[n],wd_resolution=conv_resolution[n],ws_avg=True) - gauss.reinitialize(wd_resolution=conv_resolution[n],ws_avg=True) - - X = np.copy(X0) - Y = np.copy(Y0) - - # Move select turbine around domain and compute AEP - for i in range(ny): - print('y ' + str(i+1) + ' of ' + str(ny)) - for j in range(nx): - X[idx] = xx[i,j] + 1e-8 - Y[idx] = yy[i,j] + 1e-8 - model.reinitialize(layout_x=X,layout_y=Y) - gauss.reinitialize(layout_x=X,layout_y=Y) - - aep = model.compare_aep(timer=False, display=False) - aep_flowers[n,i,j] = aep[0] - aep_park[n,i,j] = aep[1] - - aep_gauss[n,i,j] = gauss.compute_floris_aep(timer=False) - - pt = Point(xx[i,j], yy[i,j]) - - # TODO: double-check Gauss NaN issue - # if np.isnan(aep_gauss[n,i,j]): - # das - - # Mask points outside of boundary or colliding with other turbines - if np.any(np.sqrt((X[idx] - np.delete(X,idx))**2 + (Y[idx] - np.delete(Y,idx))**2) <= 126.) or not (pt.within(poly) or line.contains(pt)): - infeasible[i,j] = 1 - - # Mask and normalize by highest AEP value - infeasible = np.swapaxes(np.tile(np.expand_dims(infeasible, axis=2),len(flowers_terms)),0,2) - aep_flowers = np.ma.masked_where(infeasible,aep_flowers) - aep_park = np.ma.masked_where(infeasible,aep_park) - aep_gauss = np.ma.masked_where(infeasible,aep_gauss) - - aep_flowers /= np.expand_dims(np.amax(aep_flowers,(1,2)),(1,2)) - aep_park /= np.expand_dims(np.amax(aep_park,(1,2)),(1,2)) - aep_gauss /= np.expand_dims(np.amax(aep_gauss,(1,2)),(1,2)) - - xx /= 126. - yy /= 126. - - X0 = np.delete(X0,idx) / 126. - Y0 = np.delete(Y0,idx) / 126. - boundaries = np.array(boundaries).T / 126. - boundaries = np.append(boundaries,np.reshape(boundaries[:,0],(2,1)),axis=1) - - pickle.dump((xx, yy, aep_flowers, aep_park, aep_gauss, flowers_terms, conv_resolution, wind_rose, X0, Y0, boundaries), open(file,'wb')) diff --git a/aep_debug.py b/aep_debug.py deleted file mode 100644 index d36fa08..0000000 --- a/aep_debug.py +++ /dev/null @@ -1,11 +0,0 @@ -import model_interface as inter -import tools as tl -import numpy as np -import matplotlib.pyplot as plt - -wr = tl.load_wind_rose(7) -layout_x = 126. * np.array([0.]) -layout_y = 126. * np.array([0.]) -model = inter.AEPInterface(wr, layout_x, layout_y) -model.reinitialize(ws_avg=True) -model.compare_aep() diff --git a/aep_figures.py b/aep_figures.py deleted file mode 100644 index 96e9608..0000000 --- a/aep_figures.py +++ /dev/null @@ -1,459 +0,0 @@ -import numpy as np -import tools as tl -import matplotlib.cm as cm -import matplotlib.pyplot as plt -import matplotlib.colors as co -import pickle -import visualization as vis -from scipy.optimize import curve_fit -import matplotlib.animation as animation -import matplotlib.collections as coll -from mpl_toolkits.axes_grid1 import make_axes_locatable -from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, mark_inset - -save = False -dpi = 500 - -########################################################################### -# RANDOMIZED LAYOUT / WIND ROSE -########################################################################### -var, aep_flowers, aep_park, time_flowers, time_park = pickle.load(open('solutions/aep0.p','rb')) -wr = var[0] -nt = var[1] - -aep_flowers /= 1e9 -aep_park /= 1e9 - -# Remove outlier -idx = np.argmax(time_park) -wr = np.delete(wr,idx) -nt = np.delete(nt,idx) -aep_flowers = np.delete(aep_flowers,idx) -aep_park = np.delete(aep_park,idx) -time_flowers = np.delete(time_flowers,idx) -time_park = np.delete(time_park,idx) - - -# AEP Comparison -fig, ax = plt.subplots(1,1,figsize=(7,5)) -markers = ['o','v','^','s','P','*','X','D','p'] -cmap2 = cm.get_cmap('hot') -for i in range(9): - p = np.ma.polyfit(np.ma.masked_where(wr!=i+1,aep_park),np.ma.masked_where(wr!=i+1,aep_flowers),1) - xrange = np.array([0,np.max(np.ma.masked_where(wr!=i+1,aep_park))]) - yrange = p[0]*xrange + p[1] - im = ax.scatter(np.ma.masked_where(wr!=i+1,aep_park),np.ma.masked_where(wr!=i+1,aep_flowers),20,alpha=0.90,marker=markers[i],zorder=15,color=cmap2(0.1*i)) - ax.plot([],[],linestyle='--',markersize=5,marker=markers[i],label=str(i+1),color=cmap2(0.1*i)) - ax.plot(xrange,yrange,'--',zorder=14,color=cmap2(0.1*i)) -xlim = ax.get_xlim() -cmap = cm.get_cmap('coolwarm_r') -ax.plot([0,xlim[1]],[0,xlim[1]],'k',linewidth=2, label='0%',zorder=10) -ax.fill_between([0, xlim[1]],[0, 0.95*xlim[1]],[0, 1*xlim[1]],alpha=0.4,label='[-5%,0%]',color=cmap(0.4),zorder=1) -ax.fill_between([0, xlim[1]],[0, 0.9*xlim[1]],[0, 0.95*xlim[1]],alpha=0.4,label='[-10%,-5%]',color=cmap(0.3),zorder=1) -ax.fill_between([0, xlim[1]],[0, 0.8*xlim[1]],[0, 0.9*xlim[1]],alpha=0.4,label='[-20%,-10%]',color=cmap(0.2),zorder=1) -ax.fill_between([0, xlim[1]],[0, 0.7*xlim[1]],[0, 0.8*xlim[1]],alpha=0.4,label='[-30%,-20%]',color=cmap(0.1),zorder=1) -ax.set(xlabel='Conventional-Park AEP [GWh]',ylabel='FLOWERS AEP [GWh]',aspect='equal') -handles, labels = ax.get_legend_handles_labels() -tmp = ax.legend(handles[9:],labels[9:],bbox_to_anchor=(0.215,0.1099999,0.9,0.89),title='AEP % Diff.') -ax.add_artist(tmp) -tmp2 = ax.legend(handles[0:9],labels[0:9],loc='upper left',title='Wind Rose') -tmp2.get_frame().set_facecolor('gainsboro') -axins = zoomed_inset_axes(ax, zoom=6, loc='lower right') -for i in range(9): - p = np.ma.polyfit(np.ma.masked_where(wr!=i+1,aep_park),np.ma.masked_where(wr!=i+1,aep_flowers),1) - xrange = np.array([0,np.max(np.ma.masked_where(wr!=i+1,aep_park))]) - yrange = p[0]*xrange + p[1] - im = axins.scatter(np.ma.masked_where(wr!=i+1,aep_park),np.ma.masked_where(wr!=i+1,aep_flowers),20,alpha=0.90,marker=markers[i],zorder=15,color=cmap2(0.1*i)) - axins.plot([],[],linestyle='--',markersize=5,marker=markers[i],label=str(i+1),color=cmap2(0.1*i)) - axins.plot(xrange,yrange,'--',zorder=14,color=cmap2(0.1*i)) -xlim = axins.get_xlim() -cmap = cm.get_cmap('coolwarm_r') -axins.plot([0,xlim[1]],[0,xlim[1]],'k',linewidth=2, label='0%',zorder=10) -axins.fill_between([0, xlim[1]],[0, 0.95*xlim[1]],[0, 1*xlim[1]],alpha=0.4,label='[-5%,0%]',color=cmap(0.4),zorder=1) -axins.fill_between([0, xlim[1]],[0, 0.9*xlim[1]],[0, 0.95*xlim[1]],alpha=0.4,label='[-10%,-5%]',color=cmap(0.3),zorder=1) -axins.fill_between([0, xlim[1]],[0, 0.8*xlim[1]],[0, 0.9*xlim[1]],alpha=0.4,label='[-20%,-10%]',color=cmap(0.2),zorder=1) -axins.fill_between([0, xlim[1]],[0, 0.7*xlim[1]],[0, 0.8*xlim[1]],alpha=0.4,label='[-30%,-20%]',color=cmap(0.1),zorder=1) -axins.set(xlim=[200,1200],ylim=[200,1200],aspect='equal',xticks=[],yticks=[]) -mark_inset(ax, axins, loc1=2, loc2=4, fc="none", ec="0.5", zorder=25) - -# for i in range(9): -# p = np.ma.polyfit(np.ma.masked_where(wr!=i+1,aep_park),np.ma.masked_where(wr!=i+1,aep_flowers),1) -# xrange = np.array([0,np.max(np.ma.masked_where(wr!=i+1,aep_park))]) -# yrange = p[0]*xrange + p[1] -# im = ax[1].scatter(np.ma.masked_where(wr!=i+1,aep_park),np.ma.masked_where(wr!=i+1,aep_flowers),20,alpha=0.90,marker=markers[i],zorder=15,color=cmap2(0.1*i)) -# ax[1].plot([],[],linestyle='--',markersize=5,marker=markers[i],label=str(i+1),color=cmap2(0.1*i)) -# ax[1].plot(xrange,yrange,'--',zorder=14,color=cmap2(0.1*i)) -# xlim = ax[1].get_xlim() -# cmap = cm.get_cmap('coolwarm_r') -# ax[1].plot([0,xlim[1]],[0,xlim[1]],'k',linewidth=2, label='0%',zorder=10) -# ax[1].fill_between([0, xlim[1]],[0, 0.95*xlim[1]],[0, 0.99*xlim[1]],alpha=0.4,label='[-5%,-1%]',color=cmap(0.4),zorder=1) -# ax[1].fill_between([0, xlim[1]],[0, 0.9*xlim[1]],[0, 0.95*xlim[1]],alpha=0.4,label='[-10%,-5%]',color=cmap(0.3),zorder=1) -# ax[1].fill_between([0, xlim[1]],[0, 0.8*xlim[1]],[0, 0.9*xlim[1]],alpha=0.4,label='[-20%,-10%]',color=cmap(0.2),zorder=1) -# ax[1].fill_between([0, xlim[1]],[0, 0.7*xlim[1]],[0, 0.8*xlim[1]],alpha=0.4,label='[-30%,-20%]',color=cmap(0.1),zorder=1) -# ax[1].set(xlabel='Conventional-Park AEP [GWh]',ylabel='FLOWERS AEP [GWh]',xlim=[500,2500],ylim=[500,2500],aspect='equal',xticks=[500,1000,1500,2000,2500],yticks=[500,1000,1500,2000,2500]) -fig.tight_layout() -if save: - plt.savefig('./figures/aep_comparison.png', dpi=dpi) -# fig, ax = plt.subplots(1,1) -# markers = ['o','v','^','s','P','*','X','D','p'] -# cmap2 = cm.get_cmap('hot') -# for i in range(9): -# p = np.ma.polyfit(np.ma.masked_where(wr!=i+1,aep_park),np.ma.masked_where(wr!=i+1,aep_flowers),1) -# xrange = np.array([0,np.max(np.ma.masked_where(wr!=i+1,aep_park))]) -# yrange = p[0]*xrange + p[1] -# im = ax.scatter(np.ma.masked_where(wr!=i+1,aep_park),np.ma.masked_where(wr!=i+1,aep_flowers),30,alpha=0.85,marker=markers[i],zorder=15,color=cmap2(0.1*i)) -# ax.plot([],[],linestyle='--',markersize=5,marker=markers[i],label=str(i+1),color=cmap2(0.1*i)) -# ax.plot(xrange,yrange,'--',zorder=14,color=cmap2(0.1*i)) -# xlim = ax.get_xlim() -# cmap = cm.get_cmap('coolwarm_r') -# ax.plot([0,xlim[1]],[0,xlim[1]],'k',linewidth=2, label='0%',zorder=10) -# ax.fill_between([0, xlim[1]],[0, 0.95*xlim[1]],[0, 0.99*xlim[1]],alpha=0.4,label='[-5%,-1%]',color=cmap(0.4),zorder=1) -# ax.fill_between([0, xlim[1]],[0, 0.9*xlim[1]],[0, 0.95*xlim[1]],alpha=0.4,label='[-10%,-5%]',color=cmap(0.3),zorder=1) -# ax.fill_between([0, xlim[1]],[0, 0.8*xlim[1]],[0, 0.9*xlim[1]],alpha=0.4,label='[-20%,-10%]',color=cmap(0.2),zorder=1) -# ax.fill_between([0, xlim[1]],[0, 0.7*xlim[1]],[0, 0.8*xlim[1]],alpha=0.4,label='[-30%,-20%]',color=cmap(0.1),zorder=1) -# ax.set(xlabel='Conventional-Park AEP [GWh]',ylabel='FLOWERS AEP [GWh]') -# handles, labels = ax.get_legend_handles_labels() -# tmp = ax.legend(handles[9:],labels[9:],loc='lower right',title='AEP % Diff.') -# ax.add_artist(tmp) -# tmp = ax.legend(handles[0:9],labels[0:9],loc='upper left',title='Wind Rose') -# tmp.get_frame().set_facecolor('gainsboro') -# fig.tight_layout() -# if save: -# plt.savefig('./figures/aep_comparison.png', dpi=dpi) - -# Cost Comparison -fig, ax = plt.subplots(1,1) -im = ax.scatter(time_park/1e-3,time_flowers/1e-3,c=nt,alpha=0.7,cmap='plasma',zorder=5,vmin=2,vmax=500) -# im = ax.scatter(time_park/1e-3,time_flowers/1e-3,c=nt,cmap='plasma',zorder=5,vmin=2,vmax=500) -ax.set(xscale='log',yscale='log',aspect='equal') -plt.autoscale(False) -plt.colorbar(im,ax=ax,label='Number of Turbines',fraction=0.046,pad=0.04,shrink=0.75) -xlim = ax.get_ylim() -cmap = cm.get_cmap('Greens') -# ylim = ax.get_xlim() -ax.fill_betweenx(xlim,[xlim[0], xlim[1]],[10*xlim[0], 10*xlim[1]],alpha=0.4,label='1x-10x Speed',color=cmap(0.2)) -ax.fill_betweenx(xlim,[10*xlim[0], 10*xlim[1]],[20*xlim[0], 20*xlim[1]],alpha=0.4,label='10x-20x Speed',color=cmap(0.5)) -ax.fill_betweenx(xlim,[20*xlim[0], 20*xlim[1]],[50*xlim[0], 50*xlim[1]],alpha=0.4,label='20x-50x Speed',color=cmap(0.8)) -ax.set(xlabel='Conventional-Park Time [ms]',ylabel='FLOWERS Time [ms]') -ax.legend() -fig.tight_layout() -if save: - plt.savefig('./figures/cost_comparison.png', dpi=dpi) - -# Cost Scaling -def power_law(x, a, b): - return a*np.power(x, b) - -# time_flowers /= np.max(time_flowers) -# time_park /= np.max(time_park) - -# sorting = nt.argsort() -# idx = np.where(nt[sorting] > 30)[0][0] -# idx2 = np.where(nt[sorting] > 80)[0][0] -# p , _ = curve_fit(power_law,nt,time_flowers) -# q = np.polyfit(nt[sorting[idx:idx2]],time_park[sorting[idx:idx2]],1) -# q , _ = curve_fit(power_law,nt[sorting[idx:]],time_park[sorting[idx:]] - 0.015) -# r , _ = curve_fit(power_law,nt[sorting[:idx]],time_park[sorting[:idx]]) -num_turbs = np.linspace(2,500,endpoint=True) -num_turbs1 = np.linspace(30,80,endpoint=True) -num_turbs2 = np.linspace(80,500,endpoint=True) - -# fig, ax = plt.subplots(1,1) -# ax.scatter(nt,time_park/1e-3, s=40,linewidth=1, alpha=0.25,facecolors='none', edgecolors='tab:orange',zorder=2) -# ax.plot(num_turbs1,1.65*num_turbs1**(1)-12, '--', color='tab:orange', linewidth=2,label='Conventional: $\mathcal{O}(N^1)$') -# ax.plot(num_turbs2,0.107*num_turbs2**(1.532)+25, '--', color='tab:orange', linewidth=2,label='Conventional: $\mathcal{O}(N^1)$') -# ax.set(xlabel='Number of Terms',ylabel='AEP Evaluation Time [ms]') - -fig, ax = plt.subplots(1,1) -ax.scatter(nt,time_flowers/1e-3,s=10,linewidth=1, alpha=0.5, facecolors='tab:blue', edgecolors='tab:blue', label='FLOWERS') -ax.scatter(nt,time_park/1e-3,s=10,linewidth=1, alpha=0.5, facecolors='tab:orange', edgecolors='tab:orange', label='Conventional') -ax.plot(num_turbs, 0.8*(7.3e-4*num_turbs**(2)+0.35), '--', color='tab:blue', linewidth=3) -ax.plot(num_turbs1, 1.4*(1.65*num_turbs1**(1)-12), '--', color='tab:orange', linewidth=3) -ax.plot(num_turbs2, 0.8*(0.107*num_turbs2**(1.532)+25), '--', color='tab:orange', linewidth=3) -# ax.plot([0, 100], 1*np.array([25,1e3*np.max(time_park)]), '--', color='tab:orange', linewidth=2) -# ax.plot(range(1,101),power_law(range(1,101),p[0],p[1]),'--',linewidth=2, %(p[1])) -# ax.plot(range(31,101),power_law(range(31,101),q[0],q[1]),'--',linewidth=2,label='Conventional: $\mathcal{O}(N^{%.1f})$' %(q[1])) -# ax.plot(range(1,31),power_law(range(1,31),r[0],r[1]),'--',color="tab:orange",linewidth=2,label='') -ax.text(185,10,'$\mathcal{O}(N^2)$',color='tab:blue',fontsize=14) -ax.text(10,160,'$\mathcal{O}(N)$',color='tab:orange',fontsize=14) -ax.text(265,270,'$\mathcal{O}(N^{1.5})$',color='tab:orange',fontsize=14) -ax.set(xlabel='Number of Turbines, $N$',ylabel='AEP Evaluation Time [ms]',yscale='log') -ax.legend() -fig.tight_layout() -if save: - plt.savefig('./figures/cost_turbines.png', dpi=dpi) - -########################################################################### -# RANDOMIZED RESOLUTION -########################################################################### -var, aep_flowers, aep_park, time_flowers, time_park = pickle.load(open('solutions/aep1.p','rb')) -N_terms = var[0] -N_bins = var[1] - -# time_flowers /= np.max(time_flowers) -# time_park /= np.max(time_park) - -p = np.polyfit(time_flowers,N_terms,1) - - -# sorting = time_park.argsort() -# idx = np.where(time_park[sorting] > 0.4)[0][0] -# q = np.polyfit(time_park[sorting[idx:]],N_bins[sorting[idx:]],1) -# r = np.polyfit(time_park[sorting[:idx]],N_bins[sorting[:idx]],1) - -n_terms = np.linspace(2,180,endpoint=True) -n_bins = np.linspace(40,360,endpoint=True) - -# fig, ax = plt.subplots(1,1) -# ax.scatter(N_bins,time_park/1e-3, s=40,linewidth=1, alpha=0.25,facecolors='none', edgecolors='tab:orange',zorder=2) -# ax.plot(n_bins,0.16*n_bins**(1)+5.8e1, '--', color='tab:orange', linewidth=2,label='Conventional: $\mathcal{O}(N^1)$') -# ax.set(xlabel='Number of Terms',ylabel='AEP Evaluation Time [ms]') - -fig, ax = plt.subplots(1,1) -ax2 = ax.twinx() -ax.scatter(time_flowers/1e-3,N_terms,s=10,linewidth=1, alpha=0.5, facecolors='tab:blue', edgecolors='tab:blue',zorder=3,label='FLOWERS') -ax2.scatter(time_park/1e-3,N_bins,s=10,linewidth=1, alpha=0.5, facecolors='tab:orange', edgecolors='tab:orange',zorder=2) -ax.plot(0.7*(1.67e-1*n_terms**(1)+0.4),n_terms, '--', color='tab:blue', linewidth=3) -ax2.plot(0.7*(0.16*n_bins**(1)+5.8e1), n_bins, '--', color='tab:orange', linewidth=3) -# ax.plot(xrange,p[0]*xrange+p[1],'--',linewidth=2,color='tab:blue') -# ax2.plot(xrange2,q[0]*xrange2+q[1],'--',linewidth=2,color='tab:orange') -# ax2.plot(xrange1,r[0]*xrange1+r[1],'--',linewidth=2,color='tab:orange') -ax.set(xlabel='AEP Evaluation Time [ms]',xscale='log') -ax.set_ylabel('Number of Fourier Modes, $M$', color='tab:blue') -ax.tick_params(axis='y',labelcolor='tab:blue') -ax2.set_ylabel('Number of Wind Direction Bins, $\mathcal{D}$', color='tab:orange') -ax2.tick_params(axis='y',labelcolor='tab:orange') -ax.scatter([],[],s=10,linewidth=1, alpha=0.5, facecolors='tab:orange', edgecolors='tab:orange',label='Conventional') -ax2.text(3,107,'$\mathcal{O}(M)$',color='tab:blue',fontsize=14) -ax2.text(30,166,'$\mathcal{O}(\mathcal{D})$',color='tab:orange',fontsize=14) -ax.legend() -fig.tight_layout() -if save: - plt.savefig('./figures/cost_resolution.png', dpi=dpi) - -########################################################################### -# LAYOUT MUTATION -########################################################################### -x_all, y_all, aep_flowers, aep_park, aep_gauss, resolution, wind_rose = pickle.load(open('solutions/aep2.p','rb')) -x_all /= 126. -y_all /= 126. -N = len(x_all) - 1 -flowers_terms = resolution[0] -floris_resolution = resolution[1] - -fig = plt.figure(figsize=(9,5)) -ax0 = plt.subplot2grid((3,2),(0,0),rowspan=3) -ax1 = plt.subplot2grid((3,2),(0,1)) -ax2 = plt.subplot2grid((3,2),(1,1)) -ax3 = plt.subplot2grid((3,2),(2,1)) - -ax0.set(aspect='equal', xlim=[-3,17], ylim=[-3,17], xticks=[0,5,10,15], yticks=[0,5,10,15], xlabel='x/D', ylabel='y/D') -ax1.set(xlim=[0,N], ylim=[0.975,1.025], ylabel='Normalized AEP',xticklabels=[]) -ax2.set(xlim=[0,N], ylim=[0.975,1.025], ylabel='Normalized AEP',xticklabels=[]) -ax3.set(xlim=[0,N], ylim=[0.975,1.025], xlabel='Iteration', ylabel='Normalized AEP') - -plt.rcParams["axes.prop_cycle"] = plt.cycler("color", plt.cm.tab20c.colors) -colors = plt.rcParams['axes.prop_cycle'].by_key()['color'] -cmap = cm.get_cmap('Reds') - -patches = [] -for n in range(len(x_all[0])): - patches.append(ax0.add_patch(plt.Circle((x_all[0,n], y_all[0,n]), 1/2, color='tab:red'))) - -# line0, = ax0.plot([],[],"o",color='tab:red',markersize=12) -line1, = ax1.plot([],[],"-o",color=colors[0],markersize=3) -line2, = ax1.plot([],[],":o",color=colors[1],markersize=3) -line3, = ax1.plot([],[],"--o",color=colors[2],markersize=3) -line4, = ax2.plot([],[],"-o",color=colors[4],markersize=3) -line5, = ax2.plot([],[],":o",color=colors[5],markersize=3) -line6, = ax2.plot([],[],"--o",color=colors[6],markersize=3) -line7, = ax3.plot([],[],"-o",color=colors[8],markersize=3) -line8, = ax3.plot([],[],":o",color=colors[9],markersize=3) -line9, = ax3.plot([],[],"--o",color=colors[10],markersize=3) -divider = make_axes_locatable(ax0) -cbar_ax = divider.append_axes('top', size='5%', pad=0.05) -cbar = plt.colorbar(cm.ScalarMappable(cmap=cmap,norm=co.Normalize(vmin=0,vmax=N)),cax=cbar_ax,label='Iteration',orientation='horizontal') -cbar_ax.xaxis.set_ticks_position('top') -cbar_ax.xaxis.set_label_position('top') -fig.tight_layout() -ax1.legend(['FLOWERS: %.0f' %(flowers_terms[0]), 'FLOWERS: %.0f' %(flowers_terms[1]), 'FLOWERS: %.0f' %(flowers_terms[2])],loc='lower left',fontsize=8) -ax2.legend(['Park: $%.0f^\circ$' %(floris_resolution[0]),'Park: $%.0f^\circ$' %(floris_resolution[1]),'Park: $%.0f^\circ$' %(floris_resolution[2])],loc='lower left',fontsize=8) -ax3.legend(['Gauss: $%.0f^\circ$' %(floris_resolution[0]),'Gauss: $%.0f^\circ$' %(floris_resolution[1]),'Gauss: $%.0f^\circ$' %(floris_resolution[2])],loc='lower left',fontsize=8) - -if save: - plt.savefig('./figures/aep_mutation0.png', dpi=dpi) - -# Function to update turbine positions -def animate1(i): - for n in range(len(x_all[0])): - patches[n].center = x_all[i,n], y_all[i,n] - # line0.set_data(x_all[i], y_all[i]) - if i >= 1: - for n in range(len(x_all[0])): - ax0.plot([x_all[i-1,n],x_all[i,n]],[y_all[i-1,n],y_all[i,n]],"-o",color=cmap(i/N),markersize=3) - line1.set_data(range(i+1),aep_flowers[0,0:i+1]) - line2.set_data(range(i+1),aep_flowers[1,0:i+1]) - line3.set_data(range(i+1),aep_flowers[2,0:i+1]) - line4.set_data(range(i+1),aep_park[0,0:i+1]) - line5.set_data(range(i+1),aep_park[1,0:i+1]) - line6.set_data(range(i+1),aep_park[2,0:i+1]) - line7.set_data(range(i+1),aep_gauss[0,0:i+1]) - line8.set_data(range(i+1),aep_gauss[1,0:i+1]) - line9.set_data(range(i+1),aep_gauss[2,0:i+1]) - ax1.legend(['FLOWERS: %.0f' %(flowers_terms[0]), 'FLOWERS: %.0f' %(flowers_terms[1]), 'FLOWERS: %.0f' %(flowers_terms[2])],loc='lower left',fontsize=8) - ax2.legend(['Park: $%.0f^\circ$' %(floris_resolution[0]),'Park: $%.0f^\circ$' %(floris_resolution[1]),'Park: $%.0f^\circ$' %(floris_resolution[2])],loc='lower left',fontsize=8) - ax3.legend(['Gauss: $%.0f^\circ$' %(floris_resolution[0]),'Gauss: $%.0f^\circ$' %(floris_resolution[1]),'Gauss: $%.0f^\circ$' %(floris_resolution[2])],loc='lower left',fontsize=8) - return line1, line2, line3, line4, line5, line6, line7, line8, line9 - -# Animation -mov = animation.FuncAnimation(fig, animate1, frames=N+1, repeat=False) -if save: - mov.save('./figures/aep_mutation.gif', dpi=dpi, bitrate=500) - plt.savefig('./figures/aep_mutation.png', dpi=dpi) - -########################################################################### -# SOLUTION SPACE SMOOTHNESS -########################################################################### -xx, yy, flowers_aep, park_aep, gauss_aep, terms, conv_resolution, wind_rose, X0, Y0, boundary = pickle.load(open('solutions/aep3.p','rb')) - -boundary = np.append(boundary,np.reshape(boundary[:,0],(2,1)),axis=1) - -fig3, ax5 = plt.subplots(1,3,figsize=(11,4)) - -for i in [0,1,2]: - ax5[i].fill_between([-1,15],14,15,color='lightgray') - ax5[i].fill_between([-1,15],-1,0,color='lightgray') - ax5[i].fill_betweenx([-1,15],14,15,color='lightgray') - ax5[i].fill_betweenx([-1,15],-1,0,color='lightgray') - ax5[i].plot(boundary[0],boundary[1],linewidth=2,color="black",zorder=6) - circ1=[] - circ2=[] - for n in range(len(X0)): - circ1.append(plt.Circle((X0[n], Y0[n]), 1)) - circ2.append(plt.Circle((X0[n], Y0[n]), 1/2)) - coll1 = coll.PatchCollection(circ1, color='lightgray',zorder=5) - coll2 = coll.PatchCollection(circ2, color='black', zorder=8) - ax5[i].add_collection(coll1) - ax5[i].add_collection(coll2) - - -ax5[0].contour(xx,yy,flowers_aep[-2],levels=np.linspace(0.94,1.,50,endpoint=True),cmap='viridis',vmin=0.94,vmax=1.) -ax5[0].set( - title='FLOWERS: %.0f'%(terms[-2]), - xlabel='x/D', - ylabel='y/D', - aspect='equal' -) -ax5[1].contour(xx,yy,park_aep[4],levels=np.linspace(0.94,1.,50,endpoint=True),cmap='viridis',vmin=0.94,vmax=1.) -ax5[1].set( - title='Conventional-Park: %.0f'%(conv_resolution[4]) + '$^\circ$', - xlabel='x/D', - ylabel='y/D', - aspect='equal' -) -ax5[2].contour(xx,yy,gauss_aep[4],levels=np.linspace(0.94,1.,50,endpoint=True),cmap='viridis',vmin=0.94,vmax=1.) -ax5[2].set( - title='Conventional-Gauss: %.0f'%(conv_resolution[4]) + '$^\circ$', - xlabel='x/D', - ylabel='y/D', - aspect='equal' -) -fig3.tight_layout() -fig3.subplots_adjust(right=0.85) -cbar_ax = fig3.add_axes([0.88, 0.205, 0.02, 0.62]) -cbar = plt.colorbar(cm.ScalarMappable(norm=co.Normalize(vmin=0.94,vmax=1.)),cax=cbar_ax,label='Normalized AEP') -if save: - plt.savefig('./figures/aep_smooth.png', dpi=dpi) - -fig2, ax4 = plt.subplots(1,3,figsize=(11,4)) - -for i in [0,1,2]: - ax4[i].fill_between([-1,15],14,15,color='lightgray') - ax4[i].fill_between([-1,15],-1,0,color='lightgray') - ax4[i].fill_betweenx([-1,15],14,15,color='lightgray') - ax4[i].fill_betweenx([-1,15],-1,0,color='lightgray') - ax4[i].plot(boundary[0],boundary[1],linewidth=2,color="black",zorder=6) - circ1=[] - circ2=[] - for n in range(len(X0)): - circ1.append(plt.Circle((X0[n], Y0[n]), 1)) - circ2.append(plt.Circle((X0[n], Y0[n]), 1/2)) - coll1 = coll.PatchCollection(circ1, color='lightgray',zorder=5) - coll2 = coll.PatchCollection(circ2, color='black', zorder=8) - ax4[i].add_collection(coll1) - ax4[i].add_collection(coll2) - -ax4[0].contour(xx,yy,flowers_aep[-1],levels=np.linspace(0.94,1.,50,endpoint=True),cmap='viridis',vmin=0.94,vmax=1.) -ax4[0].set( - title='FLOWERS: %.0f'%(terms[-1]), - xlabel='x/D', - ylabel='y/D', - aspect='equal' -) -ax4[1].contour(xx,yy,park_aep[-1],levels=np.linspace(0.94,1.,50,endpoint=True),cmap='viridis',vmin=0.94,vmax=1.) -ax4[1].set( - title='Conventional-Park: %.0f'%(conv_resolution[-1]) + '$^\circ$', - xlabel='x/D', - ylabel='y/D', - aspect='equal' -) -ax4[2].contour(xx,yy,gauss_aep[-1],levels=np.linspace(0.94,1.,50,endpoint=True),cmap='viridis',vmin=0.94,vmax=1.) -ax4[2].set( - title='Conventional-Gauss: %.0f'%(conv_resolution[-1]) + '$^\circ$', - xlabel='x/D', - ylabel='y/D', - aspect='equal' -) -fig2.tight_layout() -fig2.subplots_adjust(right=0.85) -cbar_ax = fig2.add_axes([0.88, 0.205, 0.02, 0.62]) -cbar = plt.colorbar(cm.ScalarMappable(norm=co.Normalize(vmin=0.94,vmax=1.)),cax=cbar_ax,label='Normalized AEP') - -def animate2(i): - ax4[0].clear() - ax4[1].clear() - ax4[2].clear() - - for b in [0,1,2]: - ax4[b].fill_between([-1,15],14,15,color='lightgray') - ax4[b].fill_between([-1,15],-1,0,color='lightgray') - ax4[b].fill_betweenx([-1,15],14,15,color='lightgray') - ax4[b].fill_betweenx([-1,15],-1,0,color='lightgray') - ax4[b].plot(boundary[0],boundary[1],linewidth=2,color="black",zorder=6) - circ1=[] - circ2=[] - for n in range(len(X0)): - circ1.append(plt.Circle((X0[n], Y0[n]), 1)) - circ2.append(plt.Circle((X0[n], Y0[n]), 1/2)) - coll1 = coll.PatchCollection(circ1, color='lightgray',zorder=5) - coll2 = coll.PatchCollection(circ2, color='black', zorder=8) - ax4[b].add_collection(coll1) - ax4[b].add_collection(coll2) - - ax4[0].contour(xx,yy,flowers_aep[-1-i],levels=np.linspace(0.94,1.,50,endpoint=True),cmap='viridis',vmin=0.94,vmax=1.) - ax4[0].set( - title='FLOWERS: %.0f'%(terms[-1-i]), - xlabel='x/D', - ylabel='y/D', - aspect='equal' - ) - ax4[1].contour(xx,yy,park_aep[-1-i],levels=np.linspace(0.94,1.,50,endpoint=True),cmap='viridis',vmin=0.94,vmax=1.) - ax4[1].set( - title='Conventional-Park: %.0f'%(conv_resolution[-1-i]) + '$^\circ$', - xlabel='x/D', - ylabel='y/D', - aspect='equal' - ) - ax4[2].contour(xx,yy,gauss_aep[-1-i],levels=np.linspace(0.94,1.,50,endpoint=True),cmap='viridis',vmin=0.94,vmax=1.) - ax4[2].set( - title='Conventional-Gauss: %.0f'%(conv_resolution[-1-i]) + '$^\circ$', - xlabel='x/D', - ylabel='y/D', - aspect='equal' - ) - -mov2 = animation.FuncAnimation(fig2, animate2, interval=1000, frames=len(conv_resolution), repeat=True) -if save: - mov2.save('./figures/aep_smooth.gif', dpi=dpi, bitrate=500) - -plt.show() \ No newline at end of file diff --git a/array.sh b/array.sh deleted file mode 100644 index 81d8569..0000000 --- a/array.sh +++ /dev/null @@ -1,33 +0,0 @@ -#!/bin/bash -#SBATCH --account=windse -#SBATCH --job-name=multi -#SBATCH --time=2-00:00:00 -#SBATCH --nodes=1 -##SBATCH --ntasks-per-node=8 -##SBATCH --partition=debug -#SBATCH --output=output/multi.%a.out -#SBATCH --mail-user=michael.locascio@nrel.gov -#SBATCH --mail-type=BEGIN,END,FAIL - - -# go to starting directory -cd $HOME/flowers - -# load conda environment -source $HOME/.bashrc -env_flowers - -# get the JOB and SUBJOB ID -if [[ $SLURM_ARRAY_JOB_ID ]] ; then - export JOB_ID=$SLURM_ARRAY_JOB_ID - export SUB_ID=$SLURM_ARRAY_TASK_ID -else - export JOB_ID=$SLURM_JOB_ID - export SUB_ID=1 -fi - -# Run our job -srun --unbuffered -n 1 python ./multistart.py $SUB_ID > output/multi.$SUB_ID.log 2>&1 - -# submit as follows -# $ sbatch --array=0-100 array.sh \ No newline at end of file diff --git a/constraint_debug.py b/constraint_debug.py deleted file mode 100644 index cc0135b..0000000 --- a/constraint_debug.py +++ /dev/null @@ -1,152 +0,0 @@ -import numpy as np - -# User inputs -boundaries = [(0, 0),(10, 0),(10, 10),(0, 10)] -points = np.array([[-1,5,2,11],[-1,-4,5,11]]) -_nturbs = len(points[0]) -gradient=True - -_boundaries = np.array(boundaries).T -_nbounds = len(_boundaries[0]) - -# Compute edge information -_boundary_edge = np.roll(_boundaries,-1,axis=1) - _boundaries -_boundary_len = np.sqrt(_boundary_edge[0]**2 + _boundary_edge[1]**2) -_boundary_norm = np.array([_boundary_edge[1],-_boundary_edge[0]]) / _boundary_len -_boundary_int = (np.roll(_boundary_norm,1,axis=1) + _boundary_norm) / 2 - - -# Compute distances from turbines to boundary points -a = np.zeros((_nturbs,2,_nbounds)) -for i in range(_nturbs): - a[i] = np.expand_dims(points[:,i].T,axis=-1) - _boundaries - -# Compute projections -a_edge = np.sum(a*_boundary_edge, axis=1) / _boundary_len -a_int = np.sum(a*_boundary_norm, axis=1) -sigma = np.sign(np.sum(a*_boundary_int, axis=1)) - -# Initialize signed distance containers -C = np.zeros(_nturbs) -D = np.zeros(_nbounds) -if gradient: - Cx = np.zeros(_nturbs) - Cy = np.zeros(_nturbs) - -# Compute signed distance -for i in range(_nturbs): - for k in range(_nbounds): - if a_edge[i,k] < 0: - D[k] = np.sqrt(a[i,0,k]**2 + a[i,1,k]**2)*sigma[i,k] - elif a_edge[i,k] > _boundary_len[k]: - D[k] = np.sqrt(a[i,0,(k+1)%_nbounds]**2 + a[i,1,(k+1)%_nbounds]**2)*sigma[i,(k+1)%_nbounds] - else: - D[k] = a_int[i,k] - - # Select minimum distance - idx = np.argmin(np.abs(D)) - C[i] = D[idx] - - if gradient: - if a_edge[i,idx] < 0: - Cx[i] = (points[0,i] - _boundaries[0,idx]) / np.sqrt((_boundaries[0,idx]-points[0,i])**2 + (_boundaries[1,idx]-points[1,i])**2) - Cy[i] = (points[1,i] - _boundaries[1,idx]) / np.sqrt((_boundaries[0,idx]-points[0,i])**2 + (_boundaries[1,idx]-points[1,i])**2) - elif a_edge[i,idx] > _boundary_len[idx]: - Cx[i] = (points[0,i] - _boundaries[0,(idx+1)%_nbounds]) / np.sqrt((_boundaries[0,(idx+1)%_nbounds]-points[0,i])**2 + (_boundaries[1,(idx+1)%_nbounds]-points[1,i])**2) - Cy[i] = (points[1,i] - _boundaries[1,(idx+1)%_nbounds]) / np.sqrt((_boundaries[0,(idx+1)%_nbounds]-points[0,i])**2 + (_boundaries[1,(idx+1)%_nbounds]-points[1,i])**2) - else: - Cx[i] = -(_boundaries[1,idx] - _boundaries[1,(idx+1)%_nbounds]) / _boundary_len[idx] - Cy[i] = -(_boundaries[0,(idx+1)%_nbounds] - _boundaries[0,idx]) / _boundary_len[idx] - -# # Compute distances from turbines to boundary points -# a = np.zeros((_nturbs,2,_nbounds)) -# for i in range(_nturbs): -# a[i] = np.expand_dims(points[:,i].T,axis=-1) - _boundaries -# # Compute projections -# a_edge = np.sum(a*_boundary_edge, axis=1) / _boundary_len -# a_int = np.sum(a*_boundary_norm, axis=1) -# sigma = np.sign(np.sum(a*_boundary_int, axis=1)) - -# # Initialize signed distance containers -# C = np.zeros(_nturbs) -# D = np.zeros(_nbounds) -# if gradient: -# Cx = np.zeros(_nturbs) -# Cy = np.zeros(_nturbs) - -# # Compute signed distance -# for i in range(_nturbs): -# for k in range(_nbounds): -# if a_edge[i,k] < 0: -# D[k] = np.sqrt(a[i,0,k]**2 + a[i,1,k]**2)*sigma[i,k] -# elif a_edge[i,k] > _boundary_len[k]: -# D[k] = np.sqrt(a[i,0,(k+1)%_nbounds]**2 + a[i,1,(k+1)%_nbounds]**2)*sigma[i,k] -# else: -# D[k] = a_int[i,k] - -# # Select minimum distance -# idx = np.argmin(np.abs(D)) -# C[i] = D[idx] - -# if gradient: -# if a_edge[i,idx] < 0: -# Cx[i] = (points[0,i] - _boundaries[0,idx]) / np.sqrt((_boundaries[0,idx]-points[0,i])**2 + (_boundaries[1,idx]-points[1,i])**2) -# Cy[i] = (points[1,i] - _boundaries[1,idx]) / np.sqrt((_boundaries[0,idx]-points[0,i])**2 + (_boundaries[1,idx]-points[1,i])**2) -# elif a_edge[i,idx] > _boundary_len[idx]: -# Cx[i] = (points[0,i] - _boundaries[0,(idx+1)%_nbounds]) / np.sqrt((_boundaries[0,(idx+1)%_nbounds]-points[0,i])**2 + (_boundaries[1,(idx+1)%_nbounds]-points[1,i])**2) -# Cy[i] = (points[1,i] - _boundaries[1,(idx+1)%_nbounds]) / np.sqrt((_boundaries[0,(idx+1)%_nbounds]-points[0,i])**2 + (_boundaries[1,(idx+1)%_nbounds]-points[1,i])**2) -# else: -# Cx[i] = (_boundaries[1,idx] - _boundaries[1,(idx+1)%_nbounds]) / _boundary_len[idx] -# Cy[i] = (_boundaries[0,(idx+1)%_nbounds] - _boundaries[0,idx]) / _boundary_len[idx] - - -# ne = len(boundaries) -# nx = len(point[0]) - -# # Compute independent of points -# boundaries = np.array(boundaries).T -# e = np.roll(boundaries, -1) - boundaries -# e_mag = (np.sqrt(e[0]**2 + e[1]**2)) -# e = e / e_mag -# n = np.array([-e[1],e[0]]) -# q = np.roll(n, -1) + n - -# # Compute distances for each point -# a = np.zeros((nx,2,ne)) -# for i in range(nx): -# a[i] = np.expand_dims(point[:,i].T,axis=-1) - boundaries - -# a_tilde = np.sum(a*e, axis=1) -# a_hat = np.sum(a*n, axis=1) -# sigma = np.sign(np.sum(a*q, axis=1)) - -# C = np.zeros(nx) -# Cx = np.zeros(nx) -# Cy = np.zeros(nx) -# D = np.zeros(ne) -# for i in range(nx): -# for k in range(ne): -# if a_tilde[i,k] < 0: -# D[k] = np.sqrt(a[i,0,k]**2 + a[i,1,k]**2)*sigma[i,k] -# elif a_tilde[i,k] > e_mag[k]: -# D[k] = np.sqrt(a[i,0,(k+1)%ne]**2 + a[i,1,(k+1)%ne]**2)*sigma[i,k] -# else: -# D[k] = a_hat[i,k] -# idx = np.argmin(np.abs(D)) -# C[i] = D[idx] - -# if a_tilde[i,idx] < 0: -# Cx[i] = (point[0,i] - boundaries[0,idx]) / np.sqrt((boundaries[0,idx]-point[0,i])**2 + (boundaries[1,idx]-point[1,i])**2) -# Cy[i] = (point[1,i] - boundaries[1,idx]) / np.sqrt((boundaries[0,idx]-point[0,i])**2 + (boundaries[1,idx]-point[1,i])**2) -# elif a_tilde[i,idx] > e_mag[idx]: -# Cx[i] = (point[0,i] - boundaries[0,(idx+1)%ne]) / np.sqrt((boundaries[0,(idx+1)%ne]-point[0,i])**2 + (boundaries[1,(idx+1)%ne]-point[1,i])**2) -# Cy[i] = (point[1,i] - boundaries[1,(idx+1)%ne]) / np.sqrt((boundaries[0,(idx+1)%ne]-point[0,i])**2 + (boundaries[1,(idx+1)%ne]-point[1,i])**2) -# else: -# Cx[i] = (boundaries[1,idx] - boundaries[1,(idx+1)%ne]) / e_mag[idx] -# Cy[i] = (boundaries[0,(idx+1)%ne] - boundaries[0,idx]) / e_mag[idx] - - -print(C) -print(Cx) -print(Cy) - diff --git a/examples/compute_AEP.py b/examples/compute_AEP.py new file mode 100644 index 0000000..501b590 --- /dev/null +++ b/examples/compute_AEP.py @@ -0,0 +1,21 @@ + +import numpy as np +import pandas as pd + +from flowers import FlowersModel + + +# Create generic layout +D = 126. +layout_x = D * np.array([0.,0.,0.,7.,7.,7.,14.,14.,14.,21.,21.,21.,28.,28.,28.]) +layout_y = D * np.array([0.,7.,14.,0.,7.,14.,0.,7.,14.,0.,7.,14.,0.,7.,14.]) + +# Load in wind data +df = pd.read_csv('inputs/HKW_wind_rose.csv') + +# Setup FLOWERS model +flowers_model = FlowersModel(df, layout_x, layout_y) + +# Calculate the AEP +aep = flowers_model.calculate_aep() +print(aep) diff --git a/examples/inputs/HKW_wind_rose.csv b/examples/inputs/HKW_wind_rose.csv new file mode 100644 index 0000000..f608c6a --- /dev/null +++ b/examples/inputs/HKW_wind_rose.csv @@ -0,0 +1,1873 @@ +wd,ws,freq_val +0.0000000000,0.0000000000,0.0000228154 +0.0000000000,1.0000000000,0.0002224504 +0.0000000000,2.0000000000,0.0003650468 +0.0000000000,3.0000000000,0.0004848277 +0.0000000000,4.0000000000,0.0008726899 +0.0000000000,5.0000000000,0.0008669861 +0.0000000000,6.0000000000,0.0009753593 +0.0000000000,7.0000000000,0.0009810632 +0.0000000000,8.0000000000,0.0009354324 +0.0000000000,9.0000000000,0.0009753593 +0.0000000000,10.0000000000,0.0007300935 +0.0000000000,11.0000000000,0.0005817933 +0.0000000000,12.0000000000,0.0005190509 +0.0000000000,13.0000000000,0.0003935661 +0.0000000000,14.0000000000,0.0003080082 +0.0000000000,15.0000000000,0.0002509697 +0.0000000000,16.0000000000,0.0001882272 +0.0000000000,17.0000000000,0.0000741501 +0.0000000000,18.0000000000,0.0000171116 +0.0000000000,19.0000000000,0.0000228154 +0.0000000000,20.0000000000,0.0000000000 +0.0000000000,21.0000000000,0.0000057039 +0.0000000000,22.0000000000,0.0000000000 +0.0000000000,23.0000000000,0.0000000000 +0.0000000000,24.0000000000,0.0000000000 +0.0000000000,25.0000000000,0.0000000000 +5.0000000000,0.0000000000,0.0000171116 +5.0000000000,1.0000000000,0.0001996350 +5.0000000000,2.0000000000,0.0004449008 +5.0000000000,3.0000000000,0.0006445357 +5.0000000000,4.0000000000,0.0007186858 +5.0000000000,5.0000000000,0.0007871321 +5.0000000000,6.0000000000,0.0008783938 +5.0000000000,7.0000000000,0.0011407712 +5.0000000000,8.0000000000,0.0011008442 +5.0000000000,9.0000000000,0.0010209902 +5.0000000000,10.0000000000,0.0007072781 +5.0000000000,11.0000000000,0.0006331280 +5.0000000000,12.0000000000,0.0004848277 +5.0000000000,13.0000000000,0.0003593429 +5.0000000000,14.0000000000,0.0003194159 +5.0000000000,15.0000000000,0.0002053388 +5.0000000000,16.0000000000,0.0001026694 +5.0000000000,17.0000000000,0.0000627424 +5.0000000000,18.0000000000,0.0000228154 +5.0000000000,19.0000000000,0.0000000000 +5.0000000000,20.0000000000,0.0000057039 +5.0000000000,21.0000000000,0.0000000000 +5.0000000000,22.0000000000,0.0000000000 +5.0000000000,23.0000000000,0.0000000000 +5.0000000000,24.0000000000,0.0000000000 +5.0000000000,25.0000000000,0.0000000000 +10.0000000000,0.0000000000,0.0000285193 +10.0000000000,1.0000000000,0.0002167465 +10.0000000000,2.0000000000,0.0003479352 +10.0000000000,3.0000000000,0.0005361624 +10.0000000000,4.0000000000,0.0006559434 +10.0000000000,5.0000000000,0.0007985398 +10.0000000000,6.0000000000,0.0008898015 +10.0000000000,7.0000000000,0.0009867671 +10.0000000000,8.0000000000,0.0011008442 +10.0000000000,9.0000000000,0.0009753593 +10.0000000000,10.0000000000,0.0007814282 +10.0000000000,11.0000000000,0.0006217203 +10.0000000000,12.0000000000,0.0005076432 +10.0000000000,13.0000000000,0.0003707506 +10.0000000000,14.0000000000,0.0002395619 +10.0000000000,15.0000000000,0.0001197810 +10.0000000000,16.0000000000,0.0000969655 +10.0000000000,17.0000000000,0.0000228154 +10.0000000000,18.0000000000,0.0000171116 +10.0000000000,19.0000000000,0.0000000000 +10.0000000000,20.0000000000,0.0000171116 +10.0000000000,21.0000000000,0.0000000000 +10.0000000000,22.0000000000,0.0000000000 +10.0000000000,23.0000000000,0.0000000000 +10.0000000000,24.0000000000,0.0000000000 +10.0000000000,25.0000000000,0.0000000000 +15.0000000000,0.0000000000,0.0000342231 +15.0000000000,1.0000000000,0.0001996350 +15.0000000000,2.0000000000,0.0004391969 +15.0000000000,3.0000000000,0.0004391969 +15.0000000000,4.0000000000,0.0006958704 +15.0000000000,5.0000000000,0.0008384668 +15.0000000000,6.0000000000,0.0009468401 +15.0000000000,7.0000000000,0.0010209902 +15.0000000000,8.0000000000,0.0010209902 +15.0000000000,9.0000000000,0.0008898015 +15.0000000000,10.0000000000,0.0007186858 +15.0000000000,11.0000000000,0.0006274241 +15.0000000000,12.0000000000,0.0004677162 +15.0000000000,13.0000000000,0.0002966005 +15.0000000000,14.0000000000,0.0001483003 +15.0000000000,15.0000000000,0.0000684463 +15.0000000000,16.0000000000,0.0000285193 +15.0000000000,17.0000000000,0.0000342231 +15.0000000000,18.0000000000,0.0000171116 +15.0000000000,19.0000000000,0.0000399270 +15.0000000000,20.0000000000,0.0000057039 +15.0000000000,21.0000000000,0.0000114077 +15.0000000000,22.0000000000,0.0000000000 +15.0000000000,23.0000000000,0.0000000000 +15.0000000000,24.0000000000,0.0000000000 +15.0000000000,25.0000000000,0.0000000000 +20.0000000000,0.0000000000,0.0000399270 +20.0000000000,1.0000000000,0.0001597080 +20.0000000000,2.0000000000,0.0003536391 +20.0000000000,3.0000000000,0.0004848277 +20.0000000000,4.0000000000,0.0006787588 +20.0000000000,5.0000000000,0.0007415013 +20.0000000000,6.0000000000,0.0011008442 +20.0000000000,7.0000000000,0.0010209902 +20.0000000000,8.0000000000,0.0011065480 +20.0000000000,9.0000000000,0.0010837326 +20.0000000000,10.0000000000,0.0007529090 +20.0000000000,11.0000000000,0.0005532740 +20.0000000000,12.0000000000,0.0004449008 +20.0000000000,13.0000000000,0.0003308236 +20.0000000000,14.0000000000,0.0001996350 +20.0000000000,15.0000000000,0.0001254848 +20.0000000000,16.0000000000,0.0000969655 +20.0000000000,17.0000000000,0.0000513347 +20.0000000000,18.0000000000,0.0000171116 +20.0000000000,19.0000000000,0.0000000000 +20.0000000000,20.0000000000,0.0000000000 +20.0000000000,21.0000000000,0.0000000000 +20.0000000000,22.0000000000,0.0000000000 +20.0000000000,23.0000000000,0.0000000000 +20.0000000000,24.0000000000,0.0000000000 +20.0000000000,25.0000000000,0.0000000000 +25.0000000000,0.0000000000,0.0000627424 +25.0000000000,1.0000000000,0.0002053388 +25.0000000000,2.0000000000,0.0003422313 +25.0000000000,3.0000000000,0.0005304586 +25.0000000000,4.0000000000,0.0007243897 +25.0000000000,5.0000000000,0.0008384668 +25.0000000000,6.0000000000,0.0008099475 +25.0000000000,7.0000000000,0.0009981748 +25.0000000000,8.0000000000,0.0011350673 +25.0000000000,9.0000000000,0.0009810632 +25.0000000000,10.0000000000,0.0008099475 +25.0000000000,11.0000000000,0.0006559434 +25.0000000000,12.0000000000,0.0004049738 +25.0000000000,13.0000000000,0.0002966005 +25.0000000000,14.0000000000,0.0002053388 +25.0000000000,15.0000000000,0.0001197810 +25.0000000000,16.0000000000,0.0000969655 +25.0000000000,17.0000000000,0.0000513347 +25.0000000000,18.0000000000,0.0000171116 +25.0000000000,19.0000000000,0.0000114077 +25.0000000000,20.0000000000,0.0000000000 +25.0000000000,21.0000000000,0.0000000000 +25.0000000000,22.0000000000,0.0000000000 +25.0000000000,23.0000000000,0.0000000000 +25.0000000000,24.0000000000,0.0000000000 +25.0000000000,25.0000000000,0.0000000000 +30.0000000000,0.0000000000,0.0000399270 +30.0000000000,1.0000000000,0.0001825234 +30.0000000000,2.0000000000,0.0003821583 +30.0000000000,3.0000000000,0.0004734200 +30.0000000000,4.0000000000,0.0005076432 +30.0000000000,5.0000000000,0.0007757244 +30.0000000000,6.0000000000,0.0009183208 +30.0000000000,7.0000000000,0.0010381018 +30.0000000000,8.0000000000,0.0009981748 +30.0000000000,9.0000000000,0.0010495095 +30.0000000000,10.0000000000,0.0008612822 +30.0000000000,11.0000000000,0.0007015743 +30.0000000000,12.0000000000,0.0005418663 +30.0000000000,13.0000000000,0.0003479352 +30.0000000000,14.0000000000,0.0002281542 +30.0000000000,15.0000000000,0.0001140771 +30.0000000000,16.0000000000,0.0000855578 +30.0000000000,17.0000000000,0.0000285193 +30.0000000000,18.0000000000,0.0000171116 +30.0000000000,19.0000000000,0.0000000000 +30.0000000000,20.0000000000,0.0000000000 +30.0000000000,21.0000000000,0.0000000000 +30.0000000000,22.0000000000,0.0000000000 +30.0000000000,23.0000000000,0.0000000000 +30.0000000000,24.0000000000,0.0000000000 +30.0000000000,25.0000000000,0.0000000000 +35.0000000000,0.0000000000,0.0000342231 +35.0000000000,1.0000000000,0.0001711157 +35.0000000000,2.0000000000,0.0003365275 +35.0000000000,3.0000000000,0.0005418663 +35.0000000000,4.0000000000,0.0005418663 +35.0000000000,5.0000000000,0.0006844627 +35.0000000000,6.0000000000,0.0009183208 +35.0000000000,7.0000000000,0.0009639516 +35.0000000000,8.0000000000,0.0009639516 +35.0000000000,9.0000000000,0.0009297285 +35.0000000000,10.0000000000,0.0008384668 +35.0000000000,11.0000000000,0.0007415013 +35.0000000000,12.0000000000,0.0006103126 +35.0000000000,13.0000000000,0.0004391969 +35.0000000000,14.0000000000,0.0002794889 +35.0000000000,15.0000000000,0.0002281542 +35.0000000000,16.0000000000,0.0000741501 +35.0000000000,17.0000000000,0.0000399270 +35.0000000000,18.0000000000,0.0000114077 +35.0000000000,19.0000000000,0.0000000000 +35.0000000000,20.0000000000,0.0000000000 +35.0000000000,21.0000000000,0.0000000000 +35.0000000000,22.0000000000,0.0000000000 +35.0000000000,23.0000000000,0.0000000000 +35.0000000000,24.0000000000,0.0000000000 +35.0000000000,25.0000000000,0.0000000000 +40.0000000000,0.0000000000,0.0000285193 +40.0000000000,1.0000000000,0.0001825234 +40.0000000000,2.0000000000,0.0003479352 +40.0000000000,3.0000000000,0.0004791239 +40.0000000000,4.0000000000,0.0003992699 +40.0000000000,5.0000000000,0.0007472051 +40.0000000000,6.0000000000,0.0008669861 +40.0000000000,7.0000000000,0.0010381018 +40.0000000000,8.0000000000,0.0010381018 +40.0000000000,9.0000000000,0.0011407712 +40.0000000000,10.0000000000,0.0009297285 +40.0000000000,11.0000000000,0.0008955054 +40.0000000000,12.0000000000,0.0008384668 +40.0000000000,13.0000000000,0.0006103126 +40.0000000000,14.0000000000,0.0003992699 +40.0000000000,15.0000000000,0.0002623774 +40.0000000000,16.0000000000,0.0001197810 +40.0000000000,17.0000000000,0.0000456308 +40.0000000000,18.0000000000,0.0000570386 +40.0000000000,19.0000000000,0.0000171116 +40.0000000000,20.0000000000,0.0000057039 +40.0000000000,21.0000000000,0.0000000000 +40.0000000000,22.0000000000,0.0000000000 +40.0000000000,23.0000000000,0.0000000000 +40.0000000000,24.0000000000,0.0000000000 +40.0000000000,25.0000000000,0.0000000000 +45.0000000000,0.0000000000,0.0000228154 +45.0000000000,1.0000000000,0.0001368925 +45.0000000000,2.0000000000,0.0003080082 +45.0000000000,3.0000000000,0.0004391969 +45.0000000000,4.0000000000,0.0006616473 +45.0000000000,5.0000000000,0.0007700205 +45.0000000000,6.0000000000,0.0008726899 +45.0000000000,7.0000000000,0.0010095825 +45.0000000000,8.0000000000,0.0011122519 +45.0000000000,9.0000000000,0.0013175907 +45.0000000000,10.0000000000,0.0010552133 +45.0000000000,11.0000000000,0.0010038786 +45.0000000000,12.0000000000,0.0008498745 +45.0000000000,13.0000000000,0.0006844627 +45.0000000000,14.0000000000,0.0004563085 +45.0000000000,15.0000000000,0.0002794889 +45.0000000000,16.0000000000,0.0001026694 +45.0000000000,17.0000000000,0.0000399270 +45.0000000000,18.0000000000,0.0000513347 +45.0000000000,19.0000000000,0.0000171116 +45.0000000000,20.0000000000,0.0000114077 +45.0000000000,21.0000000000,0.0000000000 +45.0000000000,22.0000000000,0.0000057039 +45.0000000000,23.0000000000,0.0000000000 +45.0000000000,24.0000000000,0.0000000000 +45.0000000000,25.0000000000,0.0000000000 +50.0000000000,0.0000000000,0.0000342231 +50.0000000000,1.0000000000,0.0002452658 +50.0000000000,2.0000000000,0.0004449008 +50.0000000000,3.0000000000,0.0004677162 +50.0000000000,4.0000000000,0.0006388319 +50.0000000000,5.0000000000,0.0008612822 +50.0000000000,6.0000000000,0.0008612822 +50.0000000000,7.0000000000,0.0009525439 +50.0000000000,8.0000000000,0.0009297285 +50.0000000000,9.0000000000,0.0010951403 +50.0000000000,10.0000000000,0.0011350673 +50.0000000000,11.0000000000,0.0009981748 +50.0000000000,12.0000000000,0.0008783938 +50.0000000000,13.0000000000,0.0007186858 +50.0000000000,14.0000000000,0.0004905316 +50.0000000000,15.0000000000,0.0003536391 +50.0000000000,16.0000000000,0.0001597080 +50.0000000000,17.0000000000,0.0000513347 +50.0000000000,18.0000000000,0.0000342231 +50.0000000000,19.0000000000,0.0000285193 +50.0000000000,20.0000000000,0.0000000000 +50.0000000000,21.0000000000,0.0000000000 +50.0000000000,22.0000000000,0.0000000000 +50.0000000000,23.0000000000,0.0000000000 +50.0000000000,24.0000000000,0.0000000000 +50.0000000000,25.0000000000,0.0000000000 +55.0000000000,0.0000000000,0.0000513347 +55.0000000000,1.0000000000,0.0001996350 +55.0000000000,2.0000000000,0.0003194159 +55.0000000000,3.0000000000,0.0004848277 +55.0000000000,4.0000000000,0.0006160164 +55.0000000000,5.0000000000,0.0007928360 +55.0000000000,6.0000000000,0.0008042437 +55.0000000000,7.0000000000,0.0010152863 +55.0000000000,8.0000000000,0.0009354324 +55.0000000000,9.0000000000,0.0009411362 +55.0000000000,10.0000000000,0.0009468401 +55.0000000000,11.0000000000,0.0009012092 +55.0000000000,12.0000000000,0.0007072781 +55.0000000000,13.0000000000,0.0005989049 +55.0000000000,14.0000000000,0.0004620123 +55.0000000000,15.0000000000,0.0004563085 +55.0000000000,16.0000000000,0.0002794889 +55.0000000000,17.0000000000,0.0000912617 +55.0000000000,18.0000000000,0.0000741501 +55.0000000000,19.0000000000,0.0000285193 +55.0000000000,20.0000000000,0.0000171116 +55.0000000000,21.0000000000,0.0000114077 +55.0000000000,22.0000000000,0.0000000000 +55.0000000000,23.0000000000,0.0000000000 +55.0000000000,24.0000000000,0.0000000000 +55.0000000000,25.0000000000,0.0000000000 +60.0000000000,0.0000000000,0.0000342231 +60.0000000000,1.0000000000,0.0002509697 +60.0000000000,2.0000000000,0.0003878622 +60.0000000000,3.0000000000,0.0005361624 +60.0000000000,4.0000000000,0.0006103126 +60.0000000000,5.0000000000,0.0008213552 +60.0000000000,6.0000000000,0.0008612822 +60.0000000000,7.0000000000,0.0008441707 +60.0000000000,8.0000000000,0.0010894365 +60.0000000000,9.0000000000,0.0009981748 +60.0000000000,10.0000000000,0.0011635866 +60.0000000000,11.0000000000,0.0010209902 +60.0000000000,12.0000000000,0.0009126169 +60.0000000000,13.0000000000,0.0008156514 +60.0000000000,14.0000000000,0.0004905316 +60.0000000000,15.0000000000,0.0004563085 +60.0000000000,16.0000000000,0.0002737851 +60.0000000000,17.0000000000,0.0001368925 +60.0000000000,18.0000000000,0.0000570386 +60.0000000000,19.0000000000,0.0000228154 +60.0000000000,20.0000000000,0.0000000000 +60.0000000000,21.0000000000,0.0000000000 +60.0000000000,22.0000000000,0.0000000000 +60.0000000000,23.0000000000,0.0000000000 +60.0000000000,24.0000000000,0.0000000000 +60.0000000000,25.0000000000,0.0000000000 +65.0000000000,0.0000000000,0.0000228154 +65.0000000000,1.0000000000,0.0001996350 +65.0000000000,2.0000000000,0.0002966005 +65.0000000000,3.0000000000,0.0004848277 +65.0000000000,4.0000000000,0.0007072781 +65.0000000000,5.0000000000,0.0007700205 +65.0000000000,6.0000000000,0.0009924709 +65.0000000000,7.0000000000,0.0010266940 +65.0000000000,8.0000000000,0.0010951403 +65.0000000000,9.0000000000,0.0011122519 +65.0000000000,10.0000000000,0.0010723249 +65.0000000000,11.0000000000,0.0010438056 +65.0000000000,12.0000000000,0.0008099475 +65.0000000000,13.0000000000,0.0007871321 +65.0000000000,14.0000000000,0.0004848277 +65.0000000000,15.0000000000,0.0004449008 +65.0000000000,16.0000000000,0.0001825234 +65.0000000000,17.0000000000,0.0001540041 +65.0000000000,18.0000000000,0.0000228154 +65.0000000000,19.0000000000,0.0000000000 +65.0000000000,20.0000000000,0.0000000000 +65.0000000000,21.0000000000,0.0000000000 +65.0000000000,22.0000000000,0.0000000000 +65.0000000000,23.0000000000,0.0000000000 +65.0000000000,24.0000000000,0.0000000000 +65.0000000000,25.0000000000,0.0000000000 +70.0000000000,0.0000000000,0.0000342231 +70.0000000000,1.0000000000,0.0001254848 +70.0000000000,2.0000000000,0.0002966005 +70.0000000000,3.0000000000,0.0005703856 +70.0000000000,4.0000000000,0.0005817933 +70.0000000000,5.0000000000,0.0008156514 +70.0000000000,6.0000000000,0.0009012092 +70.0000000000,7.0000000000,0.0012263290 +70.0000000000,8.0000000000,0.0010894365 +70.0000000000,9.0000000000,0.0012092174 +70.0000000000,10.0000000000,0.0010038786 +70.0000000000,11.0000000000,0.0012548483 +70.0000000000,12.0000000000,0.0009639516 +70.0000000000,13.0000000000,0.0005874971 +70.0000000000,14.0000000000,0.0004677162 +70.0000000000,15.0000000000,0.0003878622 +70.0000000000,16.0000000000,0.0001939311 +70.0000000000,17.0000000000,0.0001254848 +70.0000000000,18.0000000000,0.0000285193 +70.0000000000,19.0000000000,0.0000057039 +70.0000000000,20.0000000000,0.0000000000 +70.0000000000,21.0000000000,0.0000114077 +70.0000000000,22.0000000000,0.0000000000 +70.0000000000,23.0000000000,0.0000000000 +70.0000000000,24.0000000000,0.0000000000 +70.0000000000,25.0000000000,0.0000000000 +75.0000000000,0.0000000000,0.0000456308 +75.0000000000,1.0000000000,0.0002167465 +75.0000000000,2.0000000000,0.0003308236 +75.0000000000,3.0000000000,0.0004677162 +75.0000000000,4.0000000000,0.0006901666 +75.0000000000,5.0000000000,0.0006616473 +75.0000000000,6.0000000000,0.0009240246 +75.0000000000,7.0000000000,0.0010438056 +75.0000000000,8.0000000000,0.0011293634 +75.0000000000,9.0000000000,0.0012719598 +75.0000000000,10.0000000000,0.0012491444 +75.0000000000,11.0000000000,0.0010095825 +75.0000000000,12.0000000000,0.0008783938 +75.0000000000,13.0000000000,0.0006160164 +75.0000000000,14.0000000000,0.0004620123 +75.0000000000,15.0000000000,0.0003080082 +75.0000000000,16.0000000000,0.0001368925 +75.0000000000,17.0000000000,0.0001197810 +75.0000000000,18.0000000000,0.0000513347 +75.0000000000,19.0000000000,0.0000285193 +75.0000000000,20.0000000000,0.0000114077 +75.0000000000,21.0000000000,0.0000114077 +75.0000000000,22.0000000000,0.0000000000 +75.0000000000,23.0000000000,0.0000000000 +75.0000000000,24.0000000000,0.0000000000 +75.0000000000,25.0000000000,0.0000000000 +80.0000000000,0.0000000000,0.0000570386 +80.0000000000,1.0000000000,0.0001711157 +80.0000000000,2.0000000000,0.0003251198 +80.0000000000,3.0000000000,0.0003935661 +80.0000000000,4.0000000000,0.0005989049 +80.0000000000,5.0000000000,0.0008042437 +80.0000000000,6.0000000000,0.0009924709 +80.0000000000,7.0000000000,0.0012263290 +80.0000000000,8.0000000000,0.0011806982 +80.0000000000,9.0000000000,0.0012662560 +80.0000000000,10.0000000000,0.0010894365 +80.0000000000,11.0000000000,0.0011122519 +80.0000000000,12.0000000000,0.0008498745 +80.0000000000,13.0000000000,0.0006445357 +80.0000000000,14.0000000000,0.0005190509 +80.0000000000,15.0000000000,0.0003251198 +80.0000000000,16.0000000000,0.0001768195 +80.0000000000,17.0000000000,0.0001083733 +80.0000000000,18.0000000000,0.0000855578 +80.0000000000,19.0000000000,0.0000285193 +80.0000000000,20.0000000000,0.0000114077 +80.0000000000,21.0000000000,0.0000114077 +80.0000000000,22.0000000000,0.0000000000 +80.0000000000,23.0000000000,0.0000000000 +80.0000000000,24.0000000000,0.0000000000 +80.0000000000,25.0000000000,0.0000000000 +85.0000000000,0.0000000000,0.0000399270 +85.0000000000,1.0000000000,0.0002053388 +85.0000000000,2.0000000000,0.0002794889 +85.0000000000,3.0000000000,0.0005247547 +85.0000000000,4.0000000000,0.0005989049 +85.0000000000,5.0000000000,0.0007472051 +85.0000000000,6.0000000000,0.0010495095 +85.0000000000,7.0000000000,0.0012149213 +85.0000000000,8.0000000000,0.0012491444 +85.0000000000,9.0000000000,0.0013632215 +85.0000000000,10.0000000000,0.0012605521 +85.0000000000,11.0000000000,0.0010495095 +85.0000000000,12.0000000000,0.0007928360 +85.0000000000,13.0000000000,0.0006160164 +85.0000000000,14.0000000000,0.0004905316 +85.0000000000,15.0000000000,0.0003251198 +85.0000000000,16.0000000000,0.0002053388 +85.0000000000,17.0000000000,0.0001540041 +85.0000000000,18.0000000000,0.0000684463 +85.0000000000,19.0000000000,0.0000570386 +85.0000000000,20.0000000000,0.0000171116 +85.0000000000,21.0000000000,0.0000114077 +85.0000000000,22.0000000000,0.0000000000 +85.0000000000,23.0000000000,0.0000000000 +85.0000000000,24.0000000000,0.0000000000 +85.0000000000,25.0000000000,0.0000000000 +90.0000000000,0.0000000000,0.0000399270 +90.0000000000,1.0000000000,0.0001311887 +90.0000000000,2.0000000000,0.0002680812 +90.0000000000,3.0000000000,0.0004563085 +90.0000000000,4.0000000000,0.0005304586 +90.0000000000,5.0000000000,0.0008498745 +90.0000000000,6.0000000000,0.0008955054 +90.0000000000,7.0000000000,0.0010438056 +90.0000000000,8.0000000000,0.0011635866 +90.0000000000,9.0000000000,0.0011236596 +90.0000000000,10.0000000000,0.0010495095 +90.0000000000,11.0000000000,0.0008555784 +90.0000000000,12.0000000000,0.0007243897 +90.0000000000,13.0000000000,0.0006046087 +90.0000000000,14.0000000000,0.0003878622 +90.0000000000,15.0000000000,0.0002167465 +90.0000000000,16.0000000000,0.0002110427 +90.0000000000,17.0000000000,0.0001996350 +90.0000000000,18.0000000000,0.0000741501 +90.0000000000,19.0000000000,0.0000741501 +90.0000000000,20.0000000000,0.0000513347 +90.0000000000,21.0000000000,0.0000114077 +90.0000000000,22.0000000000,0.0000000000 +90.0000000000,23.0000000000,0.0000000000 +90.0000000000,24.0000000000,0.0000000000 +90.0000000000,25.0000000000,0.0000000000 +95.0000000000,0.0000000000,0.0000342231 +95.0000000000,1.0000000000,0.0002680812 +95.0000000000,2.0000000000,0.0003593429 +95.0000000000,3.0000000000,0.0005076432 +95.0000000000,4.0000000000,0.0005475702 +95.0000000000,5.0000000000,0.0005932010 +95.0000000000,6.0000000000,0.0009297285 +95.0000000000,7.0000000000,0.0009468401 +95.0000000000,8.0000000000,0.0010495095 +95.0000000000,9.0000000000,0.0010266940 +95.0000000000,10.0000000000,0.0009240246 +95.0000000000,11.0000000000,0.0009354324 +95.0000000000,12.0000000000,0.0008042437 +95.0000000000,13.0000000000,0.0006844627 +95.0000000000,14.0000000000,0.0005247547 +95.0000000000,15.0000000000,0.0004391969 +95.0000000000,16.0000000000,0.0002966005 +95.0000000000,17.0000000000,0.0001311887 +95.0000000000,18.0000000000,0.0000741501 +95.0000000000,19.0000000000,0.0000627424 +95.0000000000,20.0000000000,0.0000342231 +95.0000000000,21.0000000000,0.0000000000 +95.0000000000,22.0000000000,0.0000114077 +95.0000000000,23.0000000000,0.0000000000 +95.0000000000,24.0000000000,0.0000000000 +95.0000000000,25.0000000000,0.0000000000 +100.0000000000,0.0000000000,0.0000342231 +100.0000000000,1.0000000000,0.0001711157 +100.0000000000,2.0000000000,0.0003764545 +100.0000000000,3.0000000000,0.0004620123 +100.0000000000,4.0000000000,0.0005932010 +100.0000000000,5.0000000000,0.0007871321 +100.0000000000,6.0000000000,0.0008099475 +100.0000000000,7.0000000000,0.0009069131 +100.0000000000,8.0000000000,0.0008555784 +100.0000000000,9.0000000000,0.0009411362 +100.0000000000,10.0000000000,0.0008612822 +100.0000000000,11.0000000000,0.0008555784 +100.0000000000,12.0000000000,0.0007757244 +100.0000000000,13.0000000000,0.0007357974 +100.0000000000,14.0000000000,0.0004734200 +100.0000000000,15.0000000000,0.0004049738 +100.0000000000,16.0000000000,0.0002053388 +100.0000000000,17.0000000000,0.0000741501 +100.0000000000,18.0000000000,0.0000513347 +100.0000000000,19.0000000000,0.0000114077 +100.0000000000,20.0000000000,0.0000057039 +100.0000000000,21.0000000000,0.0000000000 +100.0000000000,22.0000000000,0.0000000000 +100.0000000000,23.0000000000,0.0000000000 +100.0000000000,24.0000000000,0.0000000000 +100.0000000000,25.0000000000,0.0000000000 +105.0000000000,0.0000000000,0.0000114077 +105.0000000000,1.0000000000,0.0002509697 +105.0000000000,2.0000000000,0.0003536391 +105.0000000000,3.0000000000,0.0004848277 +105.0000000000,4.0000000000,0.0006616473 +105.0000000000,5.0000000000,0.0006616473 +105.0000000000,6.0000000000,0.0007928360 +105.0000000000,7.0000000000,0.0007985398 +105.0000000000,8.0000000000,0.0008042437 +105.0000000000,9.0000000000,0.0009069131 +105.0000000000,10.0000000000,0.0007757244 +105.0000000000,11.0000000000,0.0007072781 +105.0000000000,12.0000000000,0.0006331280 +105.0000000000,13.0000000000,0.0003935661 +105.0000000000,14.0000000000,0.0004049738 +105.0000000000,15.0000000000,0.0002338581 +105.0000000000,16.0000000000,0.0001996350 +105.0000000000,17.0000000000,0.0000741501 +105.0000000000,18.0000000000,0.0000228154 +105.0000000000,19.0000000000,0.0000000000 +105.0000000000,20.0000000000,0.0000000000 +105.0000000000,21.0000000000,0.0000000000 +105.0000000000,22.0000000000,0.0000057039 +105.0000000000,23.0000000000,0.0000000000 +105.0000000000,24.0000000000,0.0000000000 +105.0000000000,25.0000000000,0.0000000000 +110.0000000000,0.0000000000,0.0000342231 +110.0000000000,1.0000000000,0.0002110427 +110.0000000000,2.0000000000,0.0003479352 +110.0000000000,3.0000000000,0.0004391969 +110.0000000000,4.0000000000,0.0005475702 +110.0000000000,5.0000000000,0.0007072781 +110.0000000000,6.0000000000,0.0007015743 +110.0000000000,7.0000000000,0.0006901666 +110.0000000000,8.0000000000,0.0007928360 +110.0000000000,9.0000000000,0.0008726899 +110.0000000000,10.0000000000,0.0006958704 +110.0000000000,11.0000000000,0.0005703856 +110.0000000000,12.0000000000,0.0005019393 +110.0000000000,13.0000000000,0.0002966005 +110.0000000000,14.0000000000,0.0002908966 +110.0000000000,15.0000000000,0.0001654118 +110.0000000000,16.0000000000,0.0001483003 +110.0000000000,17.0000000000,0.0000399270 +110.0000000000,18.0000000000,0.0000456308 +110.0000000000,19.0000000000,0.0000057039 +110.0000000000,20.0000000000,0.0000057039 +110.0000000000,21.0000000000,0.0000000000 +110.0000000000,22.0000000000,0.0000000000 +110.0000000000,23.0000000000,0.0000000000 +110.0000000000,24.0000000000,0.0000000000 +110.0000000000,25.0000000000,0.0000000000 +115.0000000000,0.0000000000,0.0000285193 +115.0000000000,1.0000000000,0.0001425964 +115.0000000000,2.0000000000,0.0003422313 +115.0000000000,3.0000000000,0.0004506046 +115.0000000000,4.0000000000,0.0004905316 +115.0000000000,5.0000000000,0.0005361624 +115.0000000000,6.0000000000,0.0006502396 +115.0000000000,7.0000000000,0.0007928360 +115.0000000000,8.0000000000,0.0007643167 +115.0000000000,9.0000000000,0.0007757244 +115.0000000000,10.0000000000,0.0006445357 +115.0000000000,11.0000000000,0.0006673511 +115.0000000000,12.0000000000,0.0005646817 +115.0000000000,13.0000000000,0.0004106776 +115.0000000000,14.0000000000,0.0003479352 +115.0000000000,15.0000000000,0.0002509697 +115.0000000000,16.0000000000,0.0001825234 +115.0000000000,17.0000000000,0.0000513347 +115.0000000000,18.0000000000,0.0000228154 +115.0000000000,19.0000000000,0.0000171116 +115.0000000000,20.0000000000,0.0000000000 +115.0000000000,21.0000000000,0.0000000000 +115.0000000000,22.0000000000,0.0000000000 +115.0000000000,23.0000000000,0.0000000000 +115.0000000000,24.0000000000,0.0000000000 +115.0000000000,25.0000000000,0.0000000000 +120.0000000000,0.0000000000,0.0000171116 +120.0000000000,1.0000000000,0.0001311887 +120.0000000000,2.0000000000,0.0003422313 +120.0000000000,3.0000000000,0.0004049738 +120.0000000000,4.0000000000,0.0005589779 +120.0000000000,5.0000000000,0.0006046087 +120.0000000000,6.0000000000,0.0008099475 +120.0000000000,7.0000000000,0.0008783938 +120.0000000000,8.0000000000,0.0007985398 +120.0000000000,9.0000000000,0.0008955054 +120.0000000000,10.0000000000,0.0007700205 +120.0000000000,11.0000000000,0.0006388319 +120.0000000000,12.0000000000,0.0005475702 +120.0000000000,13.0000000000,0.0003593429 +120.0000000000,14.0000000000,0.0002851928 +120.0000000000,15.0000000000,0.0002908966 +120.0000000000,16.0000000000,0.0001939311 +120.0000000000,17.0000000000,0.0000741501 +120.0000000000,18.0000000000,0.0000399270 +120.0000000000,19.0000000000,0.0000057039 +120.0000000000,20.0000000000,0.0000057039 +120.0000000000,21.0000000000,0.0000000000 +120.0000000000,22.0000000000,0.0000000000 +120.0000000000,23.0000000000,0.0000000000 +120.0000000000,24.0000000000,0.0000000000 +120.0000000000,25.0000000000,0.0000000000 +125.0000000000,0.0000000000,0.0000399270 +125.0000000000,1.0000000000,0.0001483003 +125.0000000000,2.0000000000,0.0003023044 +125.0000000000,3.0000000000,0.0004563085 +125.0000000000,4.0000000000,0.0007129820 +125.0000000000,5.0000000000,0.0006673511 +125.0000000000,6.0000000000,0.0007586128 +125.0000000000,7.0000000000,0.0006901666 +125.0000000000,8.0000000000,0.0007757244 +125.0000000000,9.0000000000,0.0007529090 +125.0000000000,10.0000000000,0.0008213552 +125.0000000000,11.0000000000,0.0007015743 +125.0000000000,12.0000000000,0.0004734200 +125.0000000000,13.0000000000,0.0004106776 +125.0000000000,14.0000000000,0.0002680812 +125.0000000000,15.0000000000,0.0002110427 +125.0000000000,16.0000000000,0.0001026694 +125.0000000000,17.0000000000,0.0000627424 +125.0000000000,18.0000000000,0.0000342231 +125.0000000000,19.0000000000,0.0000228154 +125.0000000000,20.0000000000,0.0000057039 +125.0000000000,21.0000000000,0.0000057039 +125.0000000000,22.0000000000,0.0000000000 +125.0000000000,23.0000000000,0.0000000000 +125.0000000000,24.0000000000,0.0000000000 +125.0000000000,25.0000000000,0.0000000000 +130.0000000000,0.0000000000,0.0000285193 +130.0000000000,1.0000000000,0.0001768195 +130.0000000000,2.0000000000,0.0002509697 +130.0000000000,3.0000000000,0.0005361624 +130.0000000000,4.0000000000,0.0006388319 +130.0000000000,5.0000000000,0.0005589779 +130.0000000000,6.0000000000,0.0007814282 +130.0000000000,7.0000000000,0.0007529090 +130.0000000000,8.0000000000,0.0008726899 +130.0000000000,9.0000000000,0.0007357974 +130.0000000000,10.0000000000,0.0006274241 +130.0000000000,11.0000000000,0.0006502396 +130.0000000000,12.0000000000,0.0005019393 +130.0000000000,13.0000000000,0.0004277892 +130.0000000000,14.0000000000,0.0002966005 +130.0000000000,15.0000000000,0.0002224504 +130.0000000000,16.0000000000,0.0001483003 +130.0000000000,17.0000000000,0.0000570386 +130.0000000000,18.0000000000,0.0000114077 +130.0000000000,19.0000000000,0.0000114077 +130.0000000000,20.0000000000,0.0000114077 +130.0000000000,21.0000000000,0.0000000000 +130.0000000000,22.0000000000,0.0000000000 +130.0000000000,23.0000000000,0.0000000000 +130.0000000000,24.0000000000,0.0000000000 +130.0000000000,25.0000000000,0.0000000000 +135.0000000000,0.0000000000,0.0000684463 +135.0000000000,1.0000000000,0.0001768195 +135.0000000000,2.0000000000,0.0002908966 +135.0000000000,3.0000000000,0.0005019393 +135.0000000000,4.0000000000,0.0005760894 +135.0000000000,5.0000000000,0.0005874971 +135.0000000000,6.0000000000,0.0006673511 +135.0000000000,7.0000000000,0.0007814282 +135.0000000000,8.0000000000,0.0008042437 +135.0000000000,9.0000000000,0.0007357974 +135.0000000000,10.0000000000,0.0005989049 +135.0000000000,11.0000000000,0.0005589779 +135.0000000000,12.0000000000,0.0003479352 +135.0000000000,13.0000000000,0.0004049738 +135.0000000000,14.0000000000,0.0002623774 +135.0000000000,15.0000000000,0.0001425964 +135.0000000000,16.0000000000,0.0001368925 +135.0000000000,17.0000000000,0.0000513347 +135.0000000000,18.0000000000,0.0000285193 +135.0000000000,19.0000000000,0.0000114077 +135.0000000000,20.0000000000,0.0000000000 +135.0000000000,21.0000000000,0.0000000000 +135.0000000000,22.0000000000,0.0000000000 +135.0000000000,23.0000000000,0.0000000000 +135.0000000000,24.0000000000,0.0000000000 +135.0000000000,25.0000000000,0.0000000000 +140.0000000000,0.0000000000,0.0000342231 +140.0000000000,1.0000000000,0.0001425964 +140.0000000000,2.0000000000,0.0003137121 +140.0000000000,3.0000000000,0.0005589779 +140.0000000000,4.0000000000,0.0005760894 +140.0000000000,5.0000000000,0.0006331280 +140.0000000000,6.0000000000,0.0007129820 +140.0000000000,7.0000000000,0.0006388319 +140.0000000000,8.0000000000,0.0006616473 +140.0000000000,9.0000000000,0.0006559434 +140.0000000000,10.0000000000,0.0005932010 +140.0000000000,11.0000000000,0.0005361624 +140.0000000000,12.0000000000,0.0004848277 +140.0000000000,13.0000000000,0.0003593429 +140.0000000000,14.0000000000,0.0002623774 +140.0000000000,15.0000000000,0.0001083733 +140.0000000000,16.0000000000,0.0000855578 +140.0000000000,17.0000000000,0.0000513347 +140.0000000000,18.0000000000,0.0000342231 +140.0000000000,19.0000000000,0.0000114077 +140.0000000000,20.0000000000,0.0000114077 +140.0000000000,21.0000000000,0.0000000000 +140.0000000000,22.0000000000,0.0000000000 +140.0000000000,23.0000000000,0.0000000000 +140.0000000000,24.0000000000,0.0000000000 +140.0000000000,25.0000000000,0.0000000000 +145.0000000000,0.0000000000,0.0000342231 +145.0000000000,1.0000000000,0.0002224504 +145.0000000000,2.0000000000,0.0002566735 +145.0000000000,3.0000000000,0.0005190509 +145.0000000000,4.0000000000,0.0005589779 +145.0000000000,5.0000000000,0.0006046087 +145.0000000000,6.0000000000,0.0007586128 +145.0000000000,7.0000000000,0.0006274241 +145.0000000000,8.0000000000,0.0006901666 +145.0000000000,9.0000000000,0.0006217203 +145.0000000000,10.0000000000,0.0005932010 +145.0000000000,11.0000000000,0.0005475702 +145.0000000000,12.0000000000,0.0004848277 +145.0000000000,13.0000000000,0.0004049738 +145.0000000000,14.0000000000,0.0002281542 +145.0000000000,15.0000000000,0.0001026694 +145.0000000000,16.0000000000,0.0000912617 +145.0000000000,17.0000000000,0.0000570386 +145.0000000000,18.0000000000,0.0000285193 +145.0000000000,19.0000000000,0.0000228154 +145.0000000000,20.0000000000,0.0000114077 +145.0000000000,21.0000000000,0.0000057039 +145.0000000000,22.0000000000,0.0000000000 +145.0000000000,23.0000000000,0.0000000000 +145.0000000000,24.0000000000,0.0000000000 +145.0000000000,25.0000000000,0.0000000000 +150.0000000000,0.0000000000,0.0000285193 +150.0000000000,1.0000000000,0.0001768195 +150.0000000000,2.0000000000,0.0003536391 +150.0000000000,3.0000000000,0.0005989049 +150.0000000000,4.0000000000,0.0005361624 +150.0000000000,5.0000000000,0.0006160164 +150.0000000000,6.0000000000,0.0006844627 +150.0000000000,7.0000000000,0.0006730550 +150.0000000000,8.0000000000,0.0006388319 +150.0000000000,9.0000000000,0.0006787588 +150.0000000000,10.0000000000,0.0007072781 +150.0000000000,11.0000000000,0.0006844627 +150.0000000000,12.0000000000,0.0004734200 +150.0000000000,13.0000000000,0.0003080082 +150.0000000000,14.0000000000,0.0002623774 +150.0000000000,15.0000000000,0.0001540041 +150.0000000000,16.0000000000,0.0001654118 +150.0000000000,17.0000000000,0.0000798540 +150.0000000000,18.0000000000,0.0000513347 +150.0000000000,19.0000000000,0.0000171116 +150.0000000000,20.0000000000,0.0000228154 +150.0000000000,21.0000000000,0.0000057039 +150.0000000000,22.0000000000,0.0000057039 +150.0000000000,23.0000000000,0.0000000000 +150.0000000000,24.0000000000,0.0000000000 +150.0000000000,25.0000000000,0.0000000000 +155.0000000000,0.0000000000,0.0000342231 +155.0000000000,1.0000000000,0.0001996350 +155.0000000000,2.0000000000,0.0002966005 +155.0000000000,3.0000000000,0.0005133470 +155.0000000000,4.0000000000,0.0005646817 +155.0000000000,5.0000000000,0.0007186858 +155.0000000000,6.0000000000,0.0006445357 +155.0000000000,7.0000000000,0.0007357974 +155.0000000000,8.0000000000,0.0005932010 +155.0000000000,9.0000000000,0.0007928360 +155.0000000000,10.0000000000,0.0007243897 +155.0000000000,11.0000000000,0.0005760894 +155.0000000000,12.0000000000,0.0004506046 +155.0000000000,13.0000000000,0.0003650468 +155.0000000000,14.0000000000,0.0003251198 +155.0000000000,15.0000000000,0.0002680812 +155.0000000000,16.0000000000,0.0001140771 +155.0000000000,17.0000000000,0.0000855578 +155.0000000000,18.0000000000,0.0000513347 +155.0000000000,19.0000000000,0.0000171116 +155.0000000000,20.0000000000,0.0000342231 +155.0000000000,21.0000000000,0.0000114077 +155.0000000000,22.0000000000,0.0000114077 +155.0000000000,23.0000000000,0.0000000000 +155.0000000000,24.0000000000,0.0000057039 +155.0000000000,25.0000000000,0.0000000000 +160.0000000000,0.0000000000,0.0000399270 +160.0000000000,1.0000000000,0.0001197810 +160.0000000000,2.0000000000,0.0003650468 +160.0000000000,3.0000000000,0.0003935661 +160.0000000000,4.0000000000,0.0005475702 +160.0000000000,5.0000000000,0.0007415013 +160.0000000000,6.0000000000,0.0006217203 +160.0000000000,7.0000000000,0.0007015743 +160.0000000000,8.0000000000,0.0007529090 +160.0000000000,9.0000000000,0.0007357974 +160.0000000000,10.0000000000,0.0007186858 +160.0000000000,11.0000000000,0.0006787588 +160.0000000000,12.0000000000,0.0005646817 +160.0000000000,13.0000000000,0.0004905316 +160.0000000000,14.0000000000,0.0003764545 +160.0000000000,15.0000000000,0.0003080082 +160.0000000000,16.0000000000,0.0001825234 +160.0000000000,17.0000000000,0.0001311887 +160.0000000000,18.0000000000,0.0000969655 +160.0000000000,19.0000000000,0.0000684463 +160.0000000000,20.0000000000,0.0000171116 +160.0000000000,21.0000000000,0.0000171116 +160.0000000000,22.0000000000,0.0000000000 +160.0000000000,23.0000000000,0.0000000000 +160.0000000000,24.0000000000,0.0000057039 +160.0000000000,25.0000000000,0.0000000000 +165.0000000000,0.0000000000,0.0000399270 +165.0000000000,1.0000000000,0.0001825234 +165.0000000000,2.0000000000,0.0003878622 +165.0000000000,3.0000000000,0.0005361624 +165.0000000000,4.0000000000,0.0005190509 +165.0000000000,5.0000000000,0.0006103126 +165.0000000000,6.0000000000,0.0007243897 +165.0000000000,7.0000000000,0.0006331280 +165.0000000000,8.0000000000,0.0007586128 +165.0000000000,9.0000000000,0.0007472051 +165.0000000000,10.0000000000,0.0007300935 +165.0000000000,11.0000000000,0.0007871321 +165.0000000000,12.0000000000,0.0006502396 +165.0000000000,13.0000000000,0.0005703856 +165.0000000000,14.0000000000,0.0004962355 +165.0000000000,15.0000000000,0.0003194159 +165.0000000000,16.0000000000,0.0003080082 +165.0000000000,17.0000000000,0.0002053388 +165.0000000000,18.0000000000,0.0001425964 +165.0000000000,19.0000000000,0.0001197810 +165.0000000000,20.0000000000,0.0000513347 +165.0000000000,21.0000000000,0.0000342231 +165.0000000000,22.0000000000,0.0000057039 +165.0000000000,23.0000000000,0.0000000000 +165.0000000000,24.0000000000,0.0000000000 +165.0000000000,25.0000000000,0.0000000000 +170.0000000000,0.0000000000,0.0000513347 +170.0000000000,1.0000000000,0.0001311887 +170.0000000000,2.0000000000,0.0002680812 +170.0000000000,3.0000000000,0.0004734200 +170.0000000000,4.0000000000,0.0006160164 +170.0000000000,5.0000000000,0.0005989049 +170.0000000000,6.0000000000,0.0007586128 +170.0000000000,7.0000000000,0.0009012092 +170.0000000000,8.0000000000,0.0008384668 +170.0000000000,9.0000000000,0.0009411362 +170.0000000000,10.0000000000,0.0008555784 +170.0000000000,11.0000000000,0.0009126169 +170.0000000000,12.0000000000,0.0007871321 +170.0000000000,13.0000000000,0.0006559434 +170.0000000000,14.0000000000,0.0005989049 +170.0000000000,15.0000000000,0.0004962355 +170.0000000000,16.0000000000,0.0004391969 +170.0000000000,17.0000000000,0.0003536391 +170.0000000000,18.0000000000,0.0002053388 +170.0000000000,19.0000000000,0.0001425964 +170.0000000000,20.0000000000,0.0001026694 +170.0000000000,21.0000000000,0.0000570386 +170.0000000000,22.0000000000,0.0000171116 +170.0000000000,23.0000000000,0.0000057039 +170.0000000000,24.0000000000,0.0000057039 +170.0000000000,25.0000000000,0.0000057039 +175.0000000000,0.0000000000,0.0000285193 +175.0000000000,1.0000000000,0.0001711157 +175.0000000000,2.0000000000,0.0003023044 +175.0000000000,3.0000000000,0.0004220853 +175.0000000000,4.0000000000,0.0005475702 +175.0000000000,5.0000000000,0.0005532740 +175.0000000000,6.0000000000,0.0007814282 +175.0000000000,7.0000000000,0.0008783938 +175.0000000000,8.0000000000,0.0008840977 +175.0000000000,9.0000000000,0.0008840977 +175.0000000000,10.0000000000,0.0010495095 +175.0000000000,11.0000000000,0.0008840977 +175.0000000000,12.0000000000,0.0007472051 +175.0000000000,13.0000000000,0.0006046087 +175.0000000000,14.0000000000,0.0007072781 +175.0000000000,15.0000000000,0.0005475702 +175.0000000000,16.0000000000,0.0003878622 +175.0000000000,17.0000000000,0.0003650468 +175.0000000000,18.0000000000,0.0003194159 +175.0000000000,19.0000000000,0.0002053388 +175.0000000000,20.0000000000,0.0000912617 +175.0000000000,21.0000000000,0.0000912617 +175.0000000000,22.0000000000,0.0000228154 +175.0000000000,23.0000000000,0.0000114077 +175.0000000000,24.0000000000,0.0000057039 +175.0000000000,25.0000000000,0.0000171116 +180.0000000000,0.0000000000,0.0000285193 +180.0000000000,1.0000000000,0.0001597080 +180.0000000000,2.0000000000,0.0002737851 +180.0000000000,3.0000000000,0.0005247547 +180.0000000000,4.0000000000,0.0005532740 +180.0000000000,5.0000000000,0.0008612822 +180.0000000000,6.0000000000,0.0007129820 +180.0000000000,7.0000000000,0.0008156514 +180.0000000000,8.0000000000,0.0009354324 +180.0000000000,9.0000000000,0.0009411362 +180.0000000000,10.0000000000,0.0008898015 +180.0000000000,11.0000000000,0.0008213552 +180.0000000000,12.0000000000,0.0009696555 +180.0000000000,13.0000000000,0.0007700205 +180.0000000000,14.0000000000,0.0006787588 +180.0000000000,15.0000000000,0.0006160164 +180.0000000000,16.0000000000,0.0004620123 +180.0000000000,17.0000000000,0.0004220853 +180.0000000000,18.0000000000,0.0001654118 +180.0000000000,19.0000000000,0.0001939311 +180.0000000000,20.0000000000,0.0001768195 +180.0000000000,21.0000000000,0.0001311887 +180.0000000000,22.0000000000,0.0000456308 +180.0000000000,23.0000000000,0.0000171116 +180.0000000000,24.0000000000,0.0000114077 +180.0000000000,25.0000000000,0.0000228154 +185.0000000000,0.0000000000,0.0000228154 +185.0000000000,1.0000000000,0.0001540041 +185.0000000000,2.0000000000,0.0002737851 +185.0000000000,3.0000000000,0.0004449008 +185.0000000000,4.0000000000,0.0004449008 +185.0000000000,5.0000000000,0.0009126169 +185.0000000000,6.0000000000,0.0008213552 +185.0000000000,7.0000000000,0.0011008442 +185.0000000000,8.0000000000,0.0010666210 +185.0000000000,9.0000000000,0.0010951403 +185.0000000000,10.0000000000,0.0009696555 +185.0000000000,11.0000000000,0.0009411362 +185.0000000000,12.0000000000,0.0009582478 +185.0000000000,13.0000000000,0.0008156514 +185.0000000000,14.0000000000,0.0007814282 +185.0000000000,15.0000000000,0.0007814282 +185.0000000000,16.0000000000,0.0006445357 +185.0000000000,17.0000000000,0.0004677162 +185.0000000000,18.0000000000,0.0003707506 +185.0000000000,19.0000000000,0.0003650468 +185.0000000000,20.0000000000,0.0002566735 +185.0000000000,21.0000000000,0.0001825234 +185.0000000000,22.0000000000,0.0000741501 +185.0000000000,23.0000000000,0.0000342231 +185.0000000000,24.0000000000,0.0000171116 +185.0000000000,25.0000000000,0.0000171116 +190.0000000000,0.0000000000,0.0000456308 +190.0000000000,1.0000000000,0.0001425964 +190.0000000000,2.0000000000,0.0003080082 +190.0000000000,3.0000000000,0.0003764545 +190.0000000000,4.0000000000,0.0004962355 +190.0000000000,5.0000000000,0.0005932010 +190.0000000000,6.0000000000,0.0009240246 +190.0000000000,7.0000000000,0.0010323979 +190.0000000000,8.0000000000,0.0010837326 +190.0000000000,9.0000000000,0.0010438056 +190.0000000000,10.0000000000,0.0011293634 +190.0000000000,11.0000000000,0.0010552133 +190.0000000000,12.0000000000,0.0009126169 +190.0000000000,13.0000000000,0.0009981748 +190.0000000000,14.0000000000,0.0009297285 +190.0000000000,15.0000000000,0.0007472051 +190.0000000000,16.0000000000,0.0007015743 +190.0000000000,17.0000000000,0.0007072781 +190.0000000000,18.0000000000,0.0005817933 +190.0000000000,19.0000000000,0.0004791239 +190.0000000000,20.0000000000,0.0002509697 +190.0000000000,21.0000000000,0.0001711157 +190.0000000000,22.0000000000,0.0001083733 +190.0000000000,23.0000000000,0.0000456308 +190.0000000000,24.0000000000,0.0000285193 +190.0000000000,25.0000000000,0.0000171116 +195.0000000000,0.0000000000,0.0000342231 +195.0000000000,1.0000000000,0.0001825234 +195.0000000000,2.0000000000,0.0003023044 +195.0000000000,3.0000000000,0.0004563085 +195.0000000000,4.0000000000,0.0005989049 +195.0000000000,5.0000000000,0.0006787588 +195.0000000000,6.0000000000,0.0008840977 +195.0000000000,7.0000000000,0.0008898015 +195.0000000000,8.0000000000,0.0010837326 +195.0000000000,9.0000000000,0.0010038786 +195.0000000000,10.0000000000,0.0011407712 +195.0000000000,11.0000000000,0.0012548483 +195.0000000000,12.0000000000,0.0010323979 +195.0000000000,13.0000000000,0.0009696555 +195.0000000000,14.0000000000,0.0009468401 +195.0000000000,15.0000000000,0.0008099475 +195.0000000000,16.0000000000,0.0008498745 +195.0000000000,17.0000000000,0.0006274241 +195.0000000000,18.0000000000,0.0004506046 +195.0000000000,19.0000000000,0.0004791239 +195.0000000000,20.0000000000,0.0002851928 +195.0000000000,21.0000000000,0.0002110427 +195.0000000000,22.0000000000,0.0001197810 +195.0000000000,23.0000000000,0.0000513347 +195.0000000000,24.0000000000,0.0000627424 +195.0000000000,25.0000000000,0.0000342231 +200.0000000000,0.0000000000,0.0000285193 +200.0000000000,1.0000000000,0.0001368925 +200.0000000000,2.0000000000,0.0003764545 +200.0000000000,3.0000000000,0.0005247547 +200.0000000000,4.0000000000,0.0006958704 +200.0000000000,5.0000000000,0.0007015743 +200.0000000000,6.0000000000,0.0008726899 +200.0000000000,7.0000000000,0.0011635866 +200.0000000000,8.0000000000,0.0011464750 +200.0000000000,9.0000000000,0.0012263290 +200.0000000000,10.0000000000,0.0011521789 +200.0000000000,11.0000000000,0.0012320329 +200.0000000000,12.0000000000,0.0012263290 +200.0000000000,13.0000000000,0.0011407712 +200.0000000000,14.0000000000,0.0012377367 +200.0000000000,15.0000000000,0.0010894365 +200.0000000000,16.0000000000,0.0010152863 +200.0000000000,17.0000000000,0.0008669861 +200.0000000000,18.0000000000,0.0007814282 +200.0000000000,19.0000000000,0.0006103126 +200.0000000000,20.0000000000,0.0004277892 +200.0000000000,21.0000000000,0.0003365275 +200.0000000000,22.0000000000,0.0002452658 +200.0000000000,23.0000000000,0.0001425964 +200.0000000000,24.0000000000,0.0000342231 +200.0000000000,25.0000000000,0.0000285193 +205.0000000000,0.0000000000,0.0000114077 +205.0000000000,1.0000000000,0.0001540041 +205.0000000000,2.0000000000,0.0003878622 +205.0000000000,3.0000000000,0.0005361624 +205.0000000000,4.0000000000,0.0006673511 +205.0000000000,5.0000000000,0.0008270591 +205.0000000000,6.0000000000,0.0009525439 +205.0000000000,7.0000000000,0.0011578827 +205.0000000000,8.0000000000,0.0013746292 +205.0000000000,9.0000000000,0.0015457449 +205.0000000000,10.0000000000,0.0013632215 +205.0000000000,11.0000000000,0.0014658909 +205.0000000000,12.0000000000,0.0012890714 +205.0000000000,13.0000000000,0.0013461100 +205.0000000000,14.0000000000,0.0011806982 +205.0000000000,15.0000000000,0.0010780287 +205.0000000000,16.0000000000,0.0009126169 +205.0000000000,17.0000000000,0.0010951403 +205.0000000000,18.0000000000,0.0008726899 +205.0000000000,19.0000000000,0.0006388319 +205.0000000000,20.0000000000,0.0005361624 +205.0000000000,21.0000000000,0.0003308236 +205.0000000000,22.0000000000,0.0001939311 +205.0000000000,23.0000000000,0.0001254848 +205.0000000000,24.0000000000,0.0000741501 +205.0000000000,25.0000000000,0.0000456308 +210.0000000000,0.0000000000,0.0000456308 +210.0000000000,1.0000000000,0.0001483003 +210.0000000000,2.0000000000,0.0003308236 +210.0000000000,3.0000000000,0.0006331280 +210.0000000000,4.0000000000,0.0006844627 +210.0000000000,5.0000000000,0.0009696555 +210.0000000000,6.0000000000,0.0013118868 +210.0000000000,7.0000000000,0.0012320329 +210.0000000000,8.0000000000,0.0013061830 +210.0000000000,9.0000000000,0.0015343372 +210.0000000000,10.0000000000,0.0017111567 +210.0000000000,11.0000000000,0.0016255989 +210.0000000000,12.0000000000,0.0016141912 +210.0000000000,13.0000000000,0.0015685603 +210.0000000000,14.0000000000,0.0014658909 +210.0000000000,15.0000000000,0.0013689254 +210.0000000000,16.0000000000,0.0012662560 +210.0000000000,17.0000000000,0.0012662560 +210.0000000000,18.0000000000,0.0009411362 +210.0000000000,19.0000000000,0.0007300935 +210.0000000000,20.0000000000,0.0006103126 +210.0000000000,21.0000000000,0.0003479352 +210.0000000000,22.0000000000,0.0002851928 +210.0000000000,23.0000000000,0.0001597080 +210.0000000000,24.0000000000,0.0001026694 +210.0000000000,25.0000000000,0.0000399270 +215.0000000000,0.0000000000,0.0000456308 +215.0000000000,1.0000000000,0.0001825234 +215.0000000000,2.0000000000,0.0004106776 +215.0000000000,3.0000000000,0.0004734200 +215.0000000000,4.0000000000,0.0008099475 +215.0000000000,5.0000000000,0.0010266940 +215.0000000000,6.0000000000,0.0011692904 +215.0000000000,7.0000000000,0.0013575177 +215.0000000000,8.0000000000,0.0016598220 +215.0000000000,9.0000000000,0.0017510837 +215.0000000000,10.0000000000,0.0021389459 +215.0000000000,11.0000000000,0.0019221994 +215.0000000000,12.0000000000,0.0020077572 +215.0000000000,13.0000000000,0.0020305727 +215.0000000000,14.0000000000,0.0016826375 +215.0000000000,15.0000000000,0.0015856719 +215.0000000000,16.0000000000,0.0015457449 +215.0000000000,17.0000000000,0.0012890714 +215.0000000000,18.0000000000,0.0010552133 +215.0000000000,19.0000000000,0.0008156514 +215.0000000000,20.0000000000,0.0004962355 +215.0000000000,21.0000000000,0.0003992699 +215.0000000000,22.0000000000,0.0002908966 +215.0000000000,23.0000000000,0.0002623774 +215.0000000000,24.0000000000,0.0001140771 +215.0000000000,25.0000000000,0.0000570386 +220.0000000000,0.0000000000,0.0000456308 +220.0000000000,1.0000000000,0.0001996350 +220.0000000000,2.0000000000,0.0003878622 +220.0000000000,3.0000000000,0.0006046087 +220.0000000000,4.0000000000,0.0010095825 +220.0000000000,5.0000000000,0.0010266940 +220.0000000000,6.0000000000,0.0012719598 +220.0000000000,7.0000000000,0.0013289984 +220.0000000000,8.0000000000,0.0017111567 +220.0000000000,9.0000000000,0.0020647958 +220.0000000000,10.0000000000,0.0021902806 +220.0000000000,11.0000000000,0.0024412503 +220.0000000000,12.0000000000,0.0024013233 +220.0000000000,13.0000000000,0.0022073922 +220.0000000000,14.0000000000,0.0019050878 +220.0000000000,15.0000000000,0.0020362765 +220.0000000000,16.0000000000,0.0017453799 +220.0000000000,17.0000000000,0.0013860370 +220.0000000000,18.0000000000,0.0011008442 +220.0000000000,19.0000000000,0.0008384668 +220.0000000000,20.0000000000,0.0007472051 +220.0000000000,21.0000000000,0.0005989049 +220.0000000000,22.0000000000,0.0003878622 +220.0000000000,23.0000000000,0.0002737851 +220.0000000000,24.0000000000,0.0001311887 +220.0000000000,25.0000000000,0.0000627424 +225.0000000000,0.0000000000,0.0000399270 +225.0000000000,1.0000000000,0.0001882272 +225.0000000000,2.0000000000,0.0003365275 +225.0000000000,3.0000000000,0.0006388319 +225.0000000000,4.0000000000,0.0008327629 +225.0000000000,5.0000000000,0.0009069131 +225.0000000000,6.0000000000,0.0015115218 +225.0000000000,7.0000000000,0.0015343372 +225.0000000000,8.0000000000,0.0019450148 +225.0000000000,9.0000000000,0.0024013233 +225.0000000000,10.0000000000,0.0024184349 +225.0000000000,11.0000000000,0.0025781428 +225.0000000000,12.0000000000,0.0027264431 +225.0000000000,13.0000000000,0.0024412503 +225.0000000000,14.0000000000,0.0023613963 +225.0000000000,15.0000000000,0.0019906457 +225.0000000000,16.0000000000,0.0019621264 +225.0000000000,17.0000000000,0.0014316678 +225.0000000000,18.0000000000,0.0011008442 +225.0000000000,19.0000000000,0.0009069131 +225.0000000000,20.0000000000,0.0007415013 +225.0000000000,21.0000000000,0.0005817933 +225.0000000000,22.0000000000,0.0003821583 +225.0000000000,23.0000000000,0.0002110427 +225.0000000000,24.0000000000,0.0001311887 +225.0000000000,25.0000000000,0.0000627424 +230.0000000000,0.0000000000,0.0000399270 +230.0000000000,1.0000000000,0.0002395619 +230.0000000000,2.0000000000,0.0003650468 +230.0000000000,3.0000000000,0.0006046087 +230.0000000000,4.0000000000,0.0009354324 +230.0000000000,5.0000000000,0.0011578827 +230.0000000000,6.0000000000,0.0012719598 +230.0000000000,7.0000000000,0.0017453799 +230.0000000000,8.0000000000,0.0018252339 +230.0000000000,9.0000000000,0.0021731691 +230.0000000000,10.0000000000,0.0026009582 +230.0000000000,11.0000000000,0.0025268081 +230.0000000000,12.0000000000,0.0023842117 +230.0000000000,13.0000000000,0.0022986539 +230.0000000000,14.0000000000,0.0020933151 +230.0000000000,15.0000000000,0.0021674652 +230.0000000000,16.0000000000,0.0017111567 +230.0000000000,17.0000000000,0.0012605521 +230.0000000000,18.0000000000,0.0011122519 +230.0000000000,19.0000000000,0.0008270591 +230.0000000000,20.0000000000,0.0007072781 +230.0000000000,21.0000000000,0.0005703856 +230.0000000000,22.0000000000,0.0003707506 +230.0000000000,23.0000000000,0.0001768195 +230.0000000000,24.0000000000,0.0001540041 +230.0000000000,25.0000000000,0.0000684463 +235.0000000000,0.0000000000,0.0000171116 +235.0000000000,1.0000000000,0.0002395619 +235.0000000000,2.0000000000,0.0004391969 +235.0000000000,3.0000000000,0.0006901666 +235.0000000000,4.0000000000,0.0007643167 +235.0000000000,5.0000000000,0.0011749943 +235.0000000000,6.0000000000,0.0014145562 +235.0000000000,7.0000000000,0.0016084873 +235.0000000000,8.0000000000,0.0019107917 +235.0000000000,9.0000000000,0.0018594570 +235.0000000000,10.0000000000,0.0020990189 +235.0000000000,11.0000000000,0.0022587269 +235.0000000000,12.0000000000,0.0020077572 +235.0000000000,13.0000000000,0.0018309377 +235.0000000000,14.0000000000,0.0019507187 +235.0000000000,15.0000000000,0.0013004791 +235.0000000000,16.0000000000,0.0013347023 +235.0000000000,17.0000000000,0.0011521789 +235.0000000000,18.0000000000,0.0009696555 +235.0000000000,19.0000000000,0.0008783938 +235.0000000000,20.0000000000,0.0005532740 +235.0000000000,21.0000000000,0.0003764545 +235.0000000000,22.0000000000,0.0002281542 +235.0000000000,23.0000000000,0.0001939311 +235.0000000000,24.0000000000,0.0000570386 +235.0000000000,25.0000000000,0.0000342231 +240.0000000000,0.0000000000,0.0000171116 +240.0000000000,1.0000000000,0.0001882272 +240.0000000000,2.0000000000,0.0003650468 +240.0000000000,3.0000000000,0.0006388319 +240.0000000000,4.0000000000,0.0009411362 +240.0000000000,5.0000000000,0.0009696555 +240.0000000000,6.0000000000,0.0012947753 +240.0000000000,7.0000000000,0.0014887064 +240.0000000000,8.0000000000,0.0016484143 +240.0000000000,9.0000000000,0.0020819074 +240.0000000000,10.0000000000,0.0021503536 +240.0000000000,11.0000000000,0.0018993840 +240.0000000000,12.0000000000,0.0020590919 +240.0000000000,13.0000000000,0.0018024184 +240.0000000000,14.0000000000,0.0017510837 +240.0000000000,15.0000000000,0.0012719598 +240.0000000000,16.0000000000,0.0011350673 +240.0000000000,17.0000000000,0.0009411362 +240.0000000000,18.0000000000,0.0009069131 +240.0000000000,19.0000000000,0.0006217203 +240.0000000000,20.0000000000,0.0004277892 +240.0000000000,21.0000000000,0.0003992699 +240.0000000000,22.0000000000,0.0001768195 +240.0000000000,23.0000000000,0.0001368925 +240.0000000000,24.0000000000,0.0000798540 +240.0000000000,25.0000000000,0.0000570386 +245.0000000000,0.0000000000,0.0000342231 +245.0000000000,1.0000000000,0.0002680812 +245.0000000000,2.0000000000,0.0004620123 +245.0000000000,3.0000000000,0.0007243897 +245.0000000000,4.0000000000,0.0008213552 +245.0000000000,5.0000000000,0.0009753593 +245.0000000000,6.0000000000,0.0011065480 +245.0000000000,7.0000000000,0.0013232945 +245.0000000000,8.0000000000,0.0015913758 +245.0000000000,9.0000000000,0.0019336071 +245.0000000000,10.0000000000,0.0018651608 +245.0000000000,11.0000000000,0.0018252339 +245.0000000000,12.0000000000,0.0016198950 +245.0000000000,13.0000000000,0.0016255989 +245.0000000000,14.0000000000,0.0014031485 +245.0000000000,15.0000000000,0.0012206251 +245.0000000000,16.0000000000,0.0010552133 +245.0000000000,17.0000000000,0.0009354324 +245.0000000000,18.0000000000,0.0009012092 +245.0000000000,19.0000000000,0.0006046087 +245.0000000000,20.0000000000,0.0004620123 +245.0000000000,21.0000000000,0.0003137121 +245.0000000000,22.0000000000,0.0001540041 +245.0000000000,23.0000000000,0.0000456308 +245.0000000000,24.0000000000,0.0000456308 +245.0000000000,25.0000000000,0.0000627424 +250.0000000000,0.0000000000,0.0000342231 +250.0000000000,1.0000000000,0.0002053388 +250.0000000000,2.0000000000,0.0004449008 +250.0000000000,3.0000000000,0.0007357974 +250.0000000000,4.0000000000,0.0008327629 +250.0000000000,5.0000000000,0.0011464750 +250.0000000000,6.0000000000,0.0012320329 +250.0000000000,7.0000000000,0.0015286334 +250.0000000000,8.0000000000,0.0016940452 +250.0000000000,9.0000000000,0.0017396760 +250.0000000000,10.0000000000,0.0018309377 +250.0000000000,11.0000000000,0.0017624914 +250.0000000000,12.0000000000,0.0014430755 +250.0000000000,13.0000000000,0.0014487794 +250.0000000000,14.0000000000,0.0012719598 +250.0000000000,15.0000000000,0.0010723249 +250.0000000000,16.0000000000,0.0009069131 +250.0000000000,17.0000000000,0.0007871321 +250.0000000000,18.0000000000,0.0006160164 +250.0000000000,19.0000000000,0.0004905316 +250.0000000000,20.0000000000,0.0004677162 +250.0000000000,21.0000000000,0.0002623774 +250.0000000000,22.0000000000,0.0002110427 +250.0000000000,23.0000000000,0.0000684463 +250.0000000000,24.0000000000,0.0000342231 +250.0000000000,25.0000000000,0.0000342231 +255.0000000000,0.0000000000,0.0000228154 +255.0000000000,1.0000000000,0.0002338581 +255.0000000000,2.0000000000,0.0004049738 +255.0000000000,3.0000000000,0.0006103126 +255.0000000000,4.0000000000,0.0008783938 +255.0000000000,5.0000000000,0.0009924709 +255.0000000000,6.0000000000,0.0012491444 +255.0000000000,7.0000000000,0.0013347023 +255.0000000000,8.0000000000,0.0015229295 +255.0000000000,9.0000000000,0.0015799681 +255.0000000000,10.0000000000,0.0016940452 +255.0000000000,11.0000000000,0.0013917408 +255.0000000000,12.0000000000,0.0015172256 +255.0000000000,13.0000000000,0.0013632215 +255.0000000000,14.0000000000,0.0009867671 +255.0000000000,15.0000000000,0.0010495095 +255.0000000000,16.0000000000,0.0009183208 +255.0000000000,17.0000000000,0.0008441707 +255.0000000000,18.0000000000,0.0005760894 +255.0000000000,19.0000000000,0.0004334930 +255.0000000000,20.0000000000,0.0003365275 +255.0000000000,21.0000000000,0.0001254848 +255.0000000000,22.0000000000,0.0001425964 +255.0000000000,23.0000000000,0.0000627424 +255.0000000000,24.0000000000,0.0000513347 +255.0000000000,25.0000000000,0.0000285193 +260.0000000000,0.0000000000,0.0000513347 +260.0000000000,1.0000000000,0.0002737851 +260.0000000000,2.0000000000,0.0004049738 +260.0000000000,3.0000000000,0.0007985398 +260.0000000000,4.0000000000,0.0008384668 +260.0000000000,5.0000000000,0.0011008442 +260.0000000000,6.0000000000,0.0012377367 +260.0000000000,7.0000000000,0.0013404061 +260.0000000000,8.0000000000,0.0012149213 +260.0000000000,9.0000000000,0.0015514488 +260.0000000000,10.0000000000,0.0013917408 +260.0000000000,11.0000000000,0.0012719598 +260.0000000000,12.0000000000,0.0013004791 +260.0000000000,13.0000000000,0.0011578827 +260.0000000000,14.0000000000,0.0012434406 +260.0000000000,15.0000000000,0.0010095825 +260.0000000000,16.0000000000,0.0006445357 +260.0000000000,17.0000000000,0.0005589779 +260.0000000000,18.0000000000,0.0004905316 +260.0000000000,19.0000000000,0.0003422313 +260.0000000000,20.0000000000,0.0002338581 +260.0000000000,21.0000000000,0.0001654118 +260.0000000000,22.0000000000,0.0001026694 +260.0000000000,23.0000000000,0.0000456308 +260.0000000000,24.0000000000,0.0000171116 +260.0000000000,25.0000000000,0.0000228154 +265.0000000000,0.0000000000,0.0000513347 +265.0000000000,1.0000000000,0.0001825234 +265.0000000000,2.0000000000,0.0003707506 +265.0000000000,3.0000000000,0.0007700205 +265.0000000000,4.0000000000,0.0008783938 +265.0000000000,5.0000000000,0.0009639516 +265.0000000000,6.0000000000,0.0012434406 +265.0000000000,7.0000000000,0.0013175907 +265.0000000000,8.0000000000,0.0013289984 +265.0000000000,9.0000000000,0.0015001141 +265.0000000000,10.0000000000,0.0012719598 +265.0000000000,11.0000000000,0.0013461100 +265.0000000000,12.0000000000,0.0012206251 +265.0000000000,13.0000000000,0.0010552133 +265.0000000000,14.0000000000,0.0009753593 +265.0000000000,15.0000000000,0.0007985398 +265.0000000000,16.0000000000,0.0006673511 +265.0000000000,17.0000000000,0.0005304586 +265.0000000000,18.0000000000,0.0004220853 +265.0000000000,19.0000000000,0.0004506046 +265.0000000000,20.0000000000,0.0002452658 +265.0000000000,21.0000000000,0.0001711157 +265.0000000000,22.0000000000,0.0001254848 +265.0000000000,23.0000000000,0.0000456308 +265.0000000000,24.0000000000,0.0000114077 +265.0000000000,25.0000000000,0.0000057039 +270.0000000000,0.0000000000,0.0000171116 +270.0000000000,1.0000000000,0.0001825234 +270.0000000000,2.0000000000,0.0004220853 +270.0000000000,3.0000000000,0.0006730550 +270.0000000000,4.0000000000,0.0007928360 +270.0000000000,5.0000000000,0.0011350673 +270.0000000000,6.0000000000,0.0011464750 +270.0000000000,7.0000000000,0.0013175907 +270.0000000000,8.0000000000,0.0012662560 +270.0000000000,9.0000000000,0.0014145562 +270.0000000000,10.0000000000,0.0014088524 +270.0000000000,11.0000000000,0.0012263290 +270.0000000000,12.0000000000,0.0010552133 +270.0000000000,13.0000000000,0.0011008442 +270.0000000000,14.0000000000,0.0009468401 +270.0000000000,15.0000000000,0.0008555784 +270.0000000000,16.0000000000,0.0006958704 +270.0000000000,17.0000000000,0.0006046087 +270.0000000000,18.0000000000,0.0004391969 +270.0000000000,19.0000000000,0.0003194159 +270.0000000000,20.0000000000,0.0002851928 +270.0000000000,21.0000000000,0.0002110427 +270.0000000000,22.0000000000,0.0001197810 +270.0000000000,23.0000000000,0.0000285193 +270.0000000000,24.0000000000,0.0000171116 +270.0000000000,25.0000000000,0.0000171116 +275.0000000000,0.0000000000,0.0000513347 +275.0000000000,1.0000000000,0.0002509697 +275.0000000000,2.0000000000,0.0005133470 +275.0000000000,3.0000000000,0.0005646817 +275.0000000000,4.0000000000,0.0008898015 +275.0000000000,5.0000000000,0.0010495095 +275.0000000000,6.0000000000,0.0012035136 +275.0000000000,7.0000000000,0.0013061830 +275.0000000000,8.0000000000,0.0014772987 +275.0000000000,9.0000000000,0.0015058179 +275.0000000000,10.0000000000,0.0015514488 +275.0000000000,11.0000000000,0.0013118868 +275.0000000000,12.0000000000,0.0010951403 +275.0000000000,13.0000000000,0.0010438056 +275.0000000000,14.0000000000,0.0007529090 +275.0000000000,15.0000000000,0.0007243897 +275.0000000000,16.0000000000,0.0006445357 +275.0000000000,17.0000000000,0.0004734200 +275.0000000000,18.0000000000,0.0004620123 +275.0000000000,19.0000000000,0.0003707506 +275.0000000000,20.0000000000,0.0001996350 +275.0000000000,21.0000000000,0.0001425964 +275.0000000000,22.0000000000,0.0000912617 +275.0000000000,23.0000000000,0.0000285193 +275.0000000000,24.0000000000,0.0000399270 +275.0000000000,25.0000000000,0.0000171116 +280.0000000000,0.0000000000,0.0000513347 +280.0000000000,1.0000000000,0.0002167465 +280.0000000000,2.0000000000,0.0003194159 +280.0000000000,3.0000000000,0.0006616473 +280.0000000000,4.0000000000,0.0007472051 +280.0000000000,5.0000000000,0.0009069131 +280.0000000000,6.0000000000,0.0010837326 +280.0000000000,7.0000000000,0.0015343372 +280.0000000000,8.0000000000,0.0014430755 +280.0000000000,9.0000000000,0.0015571526 +280.0000000000,10.0000000000,0.0013860370 +280.0000000000,11.0000000000,0.0013860370 +280.0000000000,12.0000000000,0.0011464750 +280.0000000000,13.0000000000,0.0010095825 +280.0000000000,14.0000000000,0.0007814282 +280.0000000000,15.0000000000,0.0007700205 +280.0000000000,16.0000000000,0.0005475702 +280.0000000000,17.0000000000,0.0005361624 +280.0000000000,18.0000000000,0.0004049738 +280.0000000000,19.0000000000,0.0003422313 +280.0000000000,20.0000000000,0.0002395619 +280.0000000000,21.0000000000,0.0001825234 +280.0000000000,22.0000000000,0.0000798540 +280.0000000000,23.0000000000,0.0000399270 +280.0000000000,24.0000000000,0.0000399270 +280.0000000000,25.0000000000,0.0000171116 +285.0000000000,0.0000000000,0.0000228154 +285.0000000000,1.0000000000,0.0002224504 +285.0000000000,2.0000000000,0.0003764545 +285.0000000000,3.0000000000,0.0005304586 +285.0000000000,4.0000000000,0.0007300935 +285.0000000000,5.0000000000,0.0010266940 +285.0000000000,6.0000000000,0.0012092174 +285.0000000000,7.0000000000,0.0013175907 +285.0000000000,8.0000000000,0.0013461100 +285.0000000000,9.0000000000,0.0015970796 +285.0000000000,10.0000000000,0.0012947753 +285.0000000000,11.0000000000,0.0012035136 +285.0000000000,12.0000000000,0.0011464750 +285.0000000000,13.0000000000,0.0007928360 +285.0000000000,14.0000000000,0.0008270591 +285.0000000000,15.0000000000,0.0006901666 +285.0000000000,16.0000000000,0.0005304586 +285.0000000000,17.0000000000,0.0003878622 +285.0000000000,18.0000000000,0.0003479352 +285.0000000000,19.0000000000,0.0002966005 +285.0000000000,20.0000000000,0.0002110427 +285.0000000000,21.0000000000,0.0001540041 +285.0000000000,22.0000000000,0.0000741501 +285.0000000000,23.0000000000,0.0000399270 +285.0000000000,24.0000000000,0.0000285193 +285.0000000000,25.0000000000,0.0000000000 +290.0000000000,0.0000000000,0.0000342231 +290.0000000000,1.0000000000,0.0002395619 +290.0000000000,2.0000000000,0.0004563085 +290.0000000000,3.0000000000,0.0005190509 +290.0000000000,4.0000000000,0.0007928360 +290.0000000000,5.0000000000,0.0009525439 +290.0000000000,6.0000000000,0.0010837326 +290.0000000000,7.0000000000,0.0013347023 +290.0000000000,8.0000000000,0.0014430755 +290.0000000000,9.0000000000,0.0012320329 +290.0000000000,10.0000000000,0.0013746292 +290.0000000000,11.0000000000,0.0010780287 +290.0000000000,12.0000000000,0.0009354324 +290.0000000000,13.0000000000,0.0009582478 +290.0000000000,14.0000000000,0.0006958704 +290.0000000000,15.0000000000,0.0005703856 +290.0000000000,16.0000000000,0.0003764545 +290.0000000000,17.0000000000,0.0003878622 +290.0000000000,18.0000000000,0.0002167465 +290.0000000000,19.0000000000,0.0001882272 +290.0000000000,20.0000000000,0.0001711157 +290.0000000000,21.0000000000,0.0000969655 +290.0000000000,22.0000000000,0.0000684463 +290.0000000000,23.0000000000,0.0000228154 +290.0000000000,24.0000000000,0.0000114077 +290.0000000000,25.0000000000,0.0000171116 +295.0000000000,0.0000000000,0.0000399270 +295.0000000000,1.0000000000,0.0002281542 +295.0000000000,2.0000000000,0.0004220853 +295.0000000000,3.0000000000,0.0006388319 +295.0000000000,4.0000000000,0.0007757244 +295.0000000000,5.0000000000,0.0008955054 +295.0000000000,6.0000000000,0.0012092174 +295.0000000000,7.0000000000,0.0013232945 +295.0000000000,8.0000000000,0.0011978097 +295.0000000000,9.0000000000,0.0013860370 +295.0000000000,10.0000000000,0.0013461100 +295.0000000000,11.0000000000,0.0013404061 +295.0000000000,12.0000000000,0.0010552133 +295.0000000000,13.0000000000,0.0009525439 +295.0000000000,14.0000000000,0.0006445357 +295.0000000000,15.0000000000,0.0005361624 +295.0000000000,16.0000000000,0.0002966005 +295.0000000000,17.0000000000,0.0002851928 +295.0000000000,18.0000000000,0.0001825234 +295.0000000000,19.0000000000,0.0001483003 +295.0000000000,20.0000000000,0.0000969655 +295.0000000000,21.0000000000,0.0000570386 +295.0000000000,22.0000000000,0.0000342231 +295.0000000000,23.0000000000,0.0000171116 +295.0000000000,24.0000000000,0.0000114077 +295.0000000000,25.0000000000,0.0000114077 +300.0000000000,0.0000000000,0.0000684463 +300.0000000000,1.0000000000,0.0002281542 +300.0000000000,2.0000000000,0.0004220853 +300.0000000000,3.0000000000,0.0006046087 +300.0000000000,4.0000000000,0.0008327629 +300.0000000000,5.0000000000,0.0009126169 +300.0000000000,6.0000000000,0.0011978097 +300.0000000000,7.0000000000,0.0011978097 +300.0000000000,8.0000000000,0.0013746292 +300.0000000000,9.0000000000,0.0013974447 +300.0000000000,10.0000000000,0.0012206251 +300.0000000000,11.0000000000,0.0011635866 +300.0000000000,12.0000000000,0.0009753593 +300.0000000000,13.0000000000,0.0007757244 +300.0000000000,14.0000000000,0.0005475702 +300.0000000000,15.0000000000,0.0004563085 +300.0000000000,16.0000000000,0.0003935661 +300.0000000000,17.0000000000,0.0002794889 +300.0000000000,18.0000000000,0.0001996350 +300.0000000000,19.0000000000,0.0001368925 +300.0000000000,20.0000000000,0.0000513347 +300.0000000000,21.0000000000,0.0000627424 +300.0000000000,22.0000000000,0.0000285193 +300.0000000000,23.0000000000,0.0000228154 +300.0000000000,24.0000000000,0.0000171116 +300.0000000000,25.0000000000,0.0000057039 +305.0000000000,0.0000000000,0.0000285193 +305.0000000000,1.0000000000,0.0001996350 +305.0000000000,2.0000000000,0.0004163815 +305.0000000000,3.0000000000,0.0006217203 +305.0000000000,4.0000000000,0.0007757244 +305.0000000000,5.0000000000,0.0008783938 +305.0000000000,6.0000000000,0.0011350673 +305.0000000000,7.0000000000,0.0012890714 +305.0000000000,8.0000000000,0.0012491444 +305.0000000000,9.0000000000,0.0013917408 +305.0000000000,10.0000000000,0.0012662560 +305.0000000000,11.0000000000,0.0011293634 +305.0000000000,12.0000000000,0.0008898015 +305.0000000000,13.0000000000,0.0008612822 +305.0000000000,14.0000000000,0.0006274241 +305.0000000000,15.0000000000,0.0005190509 +305.0000000000,16.0000000000,0.0002737851 +305.0000000000,17.0000000000,0.0002509697 +305.0000000000,18.0000000000,0.0001597080 +305.0000000000,19.0000000000,0.0000627424 +305.0000000000,20.0000000000,0.0001026694 +305.0000000000,21.0000000000,0.0000513347 +305.0000000000,22.0000000000,0.0000570386 +305.0000000000,23.0000000000,0.0000057039 +305.0000000000,24.0000000000,0.0000000000 +305.0000000000,25.0000000000,0.0000000000 +310.0000000000,0.0000000000,0.0000342231 +310.0000000000,1.0000000000,0.0002167465 +310.0000000000,2.0000000000,0.0004563085 +310.0000000000,3.0000000000,0.0006103126 +310.0000000000,4.0000000000,0.0008327629 +310.0000000000,5.0000000000,0.0009753593 +310.0000000000,6.0000000000,0.0010209902 +310.0000000000,7.0000000000,0.0013803331 +310.0000000000,8.0000000000,0.0012662560 +310.0000000000,9.0000000000,0.0012491444 +310.0000000000,10.0000000000,0.0012605521 +310.0000000000,11.0000000000,0.0009924709 +310.0000000000,12.0000000000,0.0009981748 +310.0000000000,13.0000000000,0.0008156514 +310.0000000000,14.0000000000,0.0006901666 +310.0000000000,15.0000000000,0.0004848277 +310.0000000000,16.0000000000,0.0004677162 +310.0000000000,17.0000000000,0.0002395619 +310.0000000000,18.0000000000,0.0002224504 +310.0000000000,19.0000000000,0.0001311887 +310.0000000000,20.0000000000,0.0000855578 +310.0000000000,21.0000000000,0.0000798540 +310.0000000000,22.0000000000,0.0000741501 +310.0000000000,23.0000000000,0.0000285193 +310.0000000000,24.0000000000,0.0000171116 +310.0000000000,25.0000000000,0.0000000000 +315.0000000000,0.0000000000,0.0000285193 +315.0000000000,1.0000000000,0.0001939311 +315.0000000000,2.0000000000,0.0003992699 +315.0000000000,3.0000000000,0.0005532740 +315.0000000000,4.0000000000,0.0008441707 +315.0000000000,5.0000000000,0.0011008442 +315.0000000000,6.0000000000,0.0012263290 +315.0000000000,7.0000000000,0.0012491444 +315.0000000000,8.0000000000,0.0012263290 +315.0000000000,9.0000000000,0.0012320329 +315.0000000000,10.0000000000,0.0011978097 +315.0000000000,11.0000000000,0.0008213552 +315.0000000000,12.0000000000,0.0008955054 +315.0000000000,13.0000000000,0.0008213552 +315.0000000000,14.0000000000,0.0007186858 +315.0000000000,15.0000000000,0.0004848277 +315.0000000000,16.0000000000,0.0003308236 +315.0000000000,17.0000000000,0.0002966005 +315.0000000000,18.0000000000,0.0002224504 +315.0000000000,19.0000000000,0.0001425964 +315.0000000000,20.0000000000,0.0001254848 +315.0000000000,21.0000000000,0.0000912617 +315.0000000000,22.0000000000,0.0000456308 +315.0000000000,23.0000000000,0.0000399270 +315.0000000000,24.0000000000,0.0000114077 +315.0000000000,25.0000000000,0.0000000000 +320.0000000000,0.0000000000,0.0000741501 +320.0000000000,1.0000000000,0.0002509697 +320.0000000000,2.0000000000,0.0003878622 +320.0000000000,3.0000000000,0.0006559434 +320.0000000000,4.0000000000,0.0007015743 +320.0000000000,5.0000000000,0.0009582478 +320.0000000000,6.0000000000,0.0009810632 +320.0000000000,7.0000000000,0.0011806982 +320.0000000000,8.0000000000,0.0010666210 +320.0000000000,9.0000000000,0.0010152863 +320.0000000000,10.0000000000,0.0012206251 +320.0000000000,11.0000000000,0.0010152863 +320.0000000000,12.0000000000,0.0008327629 +320.0000000000,13.0000000000,0.0007129820 +320.0000000000,14.0000000000,0.0005760894 +320.0000000000,15.0000000000,0.0005532740 +320.0000000000,16.0000000000,0.0003593429 +320.0000000000,17.0000000000,0.0003137121 +320.0000000000,18.0000000000,0.0001483003 +320.0000000000,19.0000000000,0.0001425964 +320.0000000000,20.0000000000,0.0000912617 +320.0000000000,21.0000000000,0.0000456308 +320.0000000000,22.0000000000,0.0000342231 +320.0000000000,23.0000000000,0.0000228154 +320.0000000000,24.0000000000,0.0000171116 +320.0000000000,25.0000000000,0.0000000000 +325.0000000000,0.0000000000,0.0000570386 +325.0000000000,1.0000000000,0.0002110427 +325.0000000000,2.0000000000,0.0004391969 +325.0000000000,3.0000000000,0.0006160164 +325.0000000000,4.0000000000,0.0008270591 +325.0000000000,5.0000000000,0.0008955054 +325.0000000000,6.0000000000,0.0011293634 +325.0000000000,7.0000000000,0.0011236596 +325.0000000000,8.0000000000,0.0011464750 +325.0000000000,9.0000000000,0.0011179557 +325.0000000000,10.0000000000,0.0010837326 +325.0000000000,11.0000000000,0.0009639516 +325.0000000000,12.0000000000,0.0009753593 +325.0000000000,13.0000000000,0.0006958704 +325.0000000000,14.0000000000,0.0006331280 +325.0000000000,15.0000000000,0.0004734200 +325.0000000000,16.0000000000,0.0003479352 +325.0000000000,17.0000000000,0.0002281542 +325.0000000000,18.0000000000,0.0002053388 +325.0000000000,19.0000000000,0.0002110427 +325.0000000000,20.0000000000,0.0001140771 +325.0000000000,21.0000000000,0.0000627424 +325.0000000000,22.0000000000,0.0000456308 +325.0000000000,23.0000000000,0.0000114077 +325.0000000000,24.0000000000,0.0000057039 +325.0000000000,25.0000000000,0.0000000000 +330.0000000000,0.0000000000,0.0000456308 +330.0000000000,1.0000000000,0.0002167465 +330.0000000000,2.0000000000,0.0004848277 +330.0000000000,3.0000000000,0.0006673511 +330.0000000000,4.0000000000,0.0007757244 +330.0000000000,5.0000000000,0.0009525439 +330.0000000000,6.0000000000,0.0010837326 +330.0000000000,7.0000000000,0.0011978097 +330.0000000000,8.0000000000,0.0009525439 +330.0000000000,9.0000000000,0.0009240246 +330.0000000000,10.0000000000,0.0008726899 +330.0000000000,11.0000000000,0.0008498745 +330.0000000000,12.0000000000,0.0008669861 +330.0000000000,13.0000000000,0.0006901666 +330.0000000000,14.0000000000,0.0005760894 +330.0000000000,15.0000000000,0.0004620123 +330.0000000000,16.0000000000,0.0003422313 +330.0000000000,17.0000000000,0.0002566735 +330.0000000000,18.0000000000,0.0002623774 +330.0000000000,19.0000000000,0.0001254848 +330.0000000000,20.0000000000,0.0001254848 +330.0000000000,21.0000000000,0.0000399270 +330.0000000000,22.0000000000,0.0000342231 +330.0000000000,23.0000000000,0.0000000000 +330.0000000000,24.0000000000,0.0000000000 +330.0000000000,25.0000000000,0.0000000000 +335.0000000000,0.0000000000,0.0000342231 +335.0000000000,1.0000000000,0.0001654118 +335.0000000000,2.0000000000,0.0005076432 +335.0000000000,3.0000000000,0.0006673511 +335.0000000000,4.0000000000,0.0007814282 +335.0000000000,5.0000000000,0.0010666210 +335.0000000000,6.0000000000,0.0009696555 +335.0000000000,7.0000000000,0.0012320329 +335.0000000000,8.0000000000,0.0011749943 +335.0000000000,9.0000000000,0.0009924709 +335.0000000000,10.0000000000,0.0010323979 +335.0000000000,11.0000000000,0.0006901666 +335.0000000000,12.0000000000,0.0007529090 +335.0000000000,13.0000000000,0.0005589779 +335.0000000000,14.0000000000,0.0004220853 +335.0000000000,15.0000000000,0.0003365275 +335.0000000000,16.0000000000,0.0002452658 +335.0000000000,17.0000000000,0.0002053388 +335.0000000000,18.0000000000,0.0001140771 +335.0000000000,19.0000000000,0.0000912617 +335.0000000000,20.0000000000,0.0000684463 +335.0000000000,21.0000000000,0.0000171116 +335.0000000000,22.0000000000,0.0000114077 +335.0000000000,23.0000000000,0.0000000000 +335.0000000000,24.0000000000,0.0000000000 +335.0000000000,25.0000000000,0.0000000000 +340.0000000000,0.0000000000,0.0000285193 +340.0000000000,1.0000000000,0.0002395619 +340.0000000000,2.0000000000,0.0003536391 +340.0000000000,3.0000000000,0.0005019393 +340.0000000000,4.0000000000,0.0007700205 +340.0000000000,5.0000000000,0.0008669861 +340.0000000000,6.0000000000,0.0008555784 +340.0000000000,7.0000000000,0.0010837326 +340.0000000000,8.0000000000,0.0010095825 +340.0000000000,9.0000000000,0.0010266940 +340.0000000000,10.0000000000,0.0010266940 +340.0000000000,11.0000000000,0.0009012092 +340.0000000000,12.0000000000,0.0008270591 +340.0000000000,13.0000000000,0.0006559434 +340.0000000000,14.0000000000,0.0004449008 +340.0000000000,15.0000000000,0.0003365275 +340.0000000000,16.0000000000,0.0002053388 +340.0000000000,17.0000000000,0.0001425964 +340.0000000000,18.0000000000,0.0000855578 +340.0000000000,19.0000000000,0.0000855578 +340.0000000000,20.0000000000,0.0000342231 +340.0000000000,21.0000000000,0.0000057039 +340.0000000000,22.0000000000,0.0000057039 +340.0000000000,23.0000000000,0.0000000000 +340.0000000000,24.0000000000,0.0000000000 +340.0000000000,25.0000000000,0.0000000000 +345.0000000000,0.0000000000,0.0000342231 +345.0000000000,1.0000000000,0.0001825234 +345.0000000000,2.0000000000,0.0003764545 +345.0000000000,3.0000000000,0.0006103126 +345.0000000000,4.0000000000,0.0007472051 +345.0000000000,5.0000000000,0.0008726899 +345.0000000000,6.0000000000,0.0008726899 +345.0000000000,7.0000000000,0.0009297285 +345.0000000000,8.0000000000,0.0010381018 +345.0000000000,9.0000000000,0.0008099475 +345.0000000000,10.0000000000,0.0009753593 +345.0000000000,11.0000000000,0.0008612822 +345.0000000000,12.0000000000,0.0005989049 +345.0000000000,13.0000000000,0.0005703856 +345.0000000000,14.0000000000,0.0003308236 +345.0000000000,15.0000000000,0.0003194159 +345.0000000000,16.0000000000,0.0002281542 +345.0000000000,17.0000000000,0.0000798540 +345.0000000000,18.0000000000,0.0000855578 +345.0000000000,19.0000000000,0.0000399270 +345.0000000000,20.0000000000,0.0000228154 +345.0000000000,21.0000000000,0.0000000000 +345.0000000000,22.0000000000,0.0000000000 +345.0000000000,23.0000000000,0.0000000000 +345.0000000000,24.0000000000,0.0000000000 +345.0000000000,25.0000000000,0.0000000000 +350.0000000000,0.0000000000,0.0000285193 +350.0000000000,1.0000000000,0.0001711157 +350.0000000000,2.0000000000,0.0004506046 +350.0000000000,3.0000000000,0.0005532740 +350.0000000000,4.0000000000,0.0006844627 +350.0000000000,5.0000000000,0.0007015743 +350.0000000000,6.0000000000,0.0008270591 +350.0000000000,7.0000000000,0.0008840977 +350.0000000000,8.0000000000,0.0009981748 +350.0000000000,9.0000000000,0.0009297285 +350.0000000000,10.0000000000,0.0009582478 +350.0000000000,11.0000000000,0.0007015743 +350.0000000000,12.0000000000,0.0005703856 +350.0000000000,13.0000000000,0.0005817933 +350.0000000000,14.0000000000,0.0002110427 +350.0000000000,15.0000000000,0.0001654118 +350.0000000000,16.0000000000,0.0001711157 +350.0000000000,17.0000000000,0.0001026694 +350.0000000000,18.0000000000,0.0000627424 +350.0000000000,19.0000000000,0.0000228154 +350.0000000000,20.0000000000,0.0000114077 +350.0000000000,21.0000000000,0.0000114077 +350.0000000000,22.0000000000,0.0000000000 +350.0000000000,23.0000000000,0.0000000000 +350.0000000000,24.0000000000,0.0000000000 +350.0000000000,25.0000000000,0.0000000000 +355.0000000000,0.0000000000,0.0000342231 +355.0000000000,1.0000000000,0.0002167465 +355.0000000000,2.0000000000,0.0003650468 +355.0000000000,3.0000000000,0.0005418663 +355.0000000000,4.0000000000,0.0006160164 +355.0000000000,5.0000000000,0.0008955054 +355.0000000000,6.0000000000,0.0008498745 +355.0000000000,7.0000000000,0.0009126169 +355.0000000000,8.0000000000,0.0010438056 +355.0000000000,9.0000000000,0.0008898015 +355.0000000000,10.0000000000,0.0009240246 +355.0000000000,11.0000000000,0.0007415013 +355.0000000000,12.0000000000,0.0005874971 +355.0000000000,13.0000000000,0.0003365275 +355.0000000000,14.0000000000,0.0003137121 +355.0000000000,15.0000000000,0.0002338581 +355.0000000000,16.0000000000,0.0001425964 +355.0000000000,17.0000000000,0.0001140771 +355.0000000000,18.0000000000,0.0000627424 +355.0000000000,19.0000000000,0.0000114077 +355.0000000000,20.0000000000,0.0000171116 +355.0000000000,21.0000000000,0.0000000000 +355.0000000000,22.0000000000,0.0000057039 +355.0000000000,23.0000000000,0.0000000000 +355.0000000000,24.0000000000,0.0000000000 +355.0000000000,25.0000000000,0.0000000000 diff --git a/flowers/__init__.py b/flowers/__init__.py new file mode 100644 index 0000000..03bb1b9 --- /dev/null +++ b/flowers/__init__.py @@ -0,0 +1 @@ +from .flowers_model import FlowersModel diff --git a/flowers_interface.py b/flowers/flowers_model.py similarity index 73% rename from flowers_interface.py rename to flowers/flowers_model.py index fc77977..1c13502 100644 --- a/flowers_interface.py +++ b/flowers/flowers_model.py @@ -1,14 +1,15 @@ -# FLOWERS -# Michael LoCascio +import copy import numpy as np import pandas as pd -import tools as tl -class FlowersInterface(): +import flowers.tools as tl + + +class FlowersModel(): """ - Flowers is a high-level user interface to the FLOWERS AEP model. + FlowersModel is a high-level user interface to the FLOWERS AEP model. Args: wind_rose (pandas.DataFrame): A dataframe for the wind rose in the FLORIS @@ -40,49 +41,53 @@ def __init__(self, wind_rose, layout_x, layout_y, num_terms=0, k=0.05, turbine=N self.turbine = 'nrel_5MW' self.D = 126. self.U = 25.0 - + if turbine == 'iea_22MW': + self.turbine = 'iea_22MW' + self.D = 283.2 + self.U = 25.0 + self._fourier_coefficients(num_terms=num_terms) - + def reinitialize(self, wind_rose=None, layout_x=None, layout_y=None, num_terms=None, k=None): if wind_rose is not None: self.wind_rose = wind_rose self._fourier_coefficients(num_terms=num_terms) - + if num_terms is not None: self._fourier_coefficients(num_terms=num_terms) - + if layout_x is not None: self.layout_x = layout_x - + if layout_y is not None: self.layout_y = layout_y - + if k is not None: self.k = k - + ########################################################################### # User functions ########################################################################### def get_layout(self): return self.layout_x, self.layout_y - + def get_wind_rose(self): return self.wind_rose - + def get_num_modes(self): return len(self.fs) - + def calculate_aep(self, gradient=False): """ Compute farm AEP (and Cartesian gradients) for the given layout and wind rose. - + Returns: aep (float): farm AEP [Wh] gradient (numpy.array(float)): (dAEP/dx, dAEP/dy) for each turbine [Wh/m] """ - + # Power component from freestream u0 = self.fs.c[0] @@ -99,16 +104,17 @@ def calculate_aep(self, gradient=False): mask_val = self.fs.c[0] # Critical polar angle of wake edge (as a function of distance from turbine) - theta_c = np.arctan( - (1 / (2*R) + self.k * np.sqrt(1 + self.k**2 - (2*R)**(-2))) - / (-self.k / (2*R) + np.sqrt(1 + self.k**2 - (2*R)**(-2))) - ) / (2 * np.pi) - theta_c = np.nan_to_num(theta_c) - + with np.errstate(divide='ignore', invalid='ignore'): + theta_c = np.arctan( + (1 / (2*R) + self.k * np.sqrt(1 + self.k**2 - (2*R)**(-2))) + / (-self.k / (2*R) + np.sqrt(1 + self.k**2 - (2*R)**(-2))) + ) / (2 * np.pi) + theta_c = np.nan_to_num(theta_c) + # Contribution from zero-frequency Fourier mode du = self.fs.a[0] * theta_c / (2 * self.k * R + 1)**2 * ( 1 + (8 * np.pi**2 * theta_c**2 * self.k * R) / (3 * (2 * self.k * R + 1))) - + # Initialize gradient and calculate zero-frequency modes if gradient == True: grad = np.zeros((len(self.layout_x),2)) @@ -122,35 +128,27 @@ def calculate_aep(self, gradient=False): 3 * self.fs.a[0] * (1 + 2 * self.k * R) * (1 + 2 * self.k * R + 8 * np.pi**2 * self.k * R * theta_c**2) * dtdr) / ( 3 * (1 + 2*self.k*R)**4) - # Reshape variables for vectorized calculations - m = np.arange(1, len(self.fs.b)) - a = np.swapaxes(np.tile(np.expand_dims(self.fs.a[1:], axis=(1,2)),np.shape(R.T)),0,2) - b = np.swapaxes(np.tile(np.expand_dims(self.fs.b[1:], axis=(1,2)),np.shape(R.T)),0,2) - R = np.tile(np.expand_dims(R, axis=2),len(m)) - THETA = np.tile(np.expand_dims(THETA, axis=2),len(m)) - theta_c = np.tile(np.expand_dims(theta_c, axis=2),len(m)) - - # Vectorized contribution of higher Fourier modes - du += np.sum((1 / (np.pi * m * (2 * self.k * R + 1)**2) * ( - a * np.cos(2 * np.pi * m * THETA) + b * np.sin(2 * np.pi * m * THETA)) * ( + for m in np.arange(1,len(self.fs.b)): + du += (1 / (np.pi * m * (2 * self.k * R + 1)**2) * ( + self.fs.a[m] * np.cos(2 * np.pi * m * THETA) + self.fs.b[m] * np.sin(2 * np.pi * m * THETA)) * ( np.sin(2 * np.pi * m * theta_c) + 2 * self.k * R / (m**2 * (2 * self.k * R + 1)) * ( - ((2 * np.pi * theta_c * m)**2 - 2) * np.sin(2 * np.pi * m * theta_c) + 4*np.pi*m*theta_c*np.cos(2 * np.pi * m * theta_c)))), axis=2) + ((2 * np.pi * theta_c * m)**2 - 2) * np.sin(2 * np.pi * m * theta_c) + 4*np.pi*m*theta_c*np.cos(2 * np.pi * m * theta_c)))) if gradient==True: - dtdr = np.tile(np.expand_dims(dtdr, axis=2),len(m)) - - # Higher Fourier modes of change in power deficit wrt angle - dpdt = np.sum((2 / (2 * self.k * R + 1)**2 * ( - b * np.cos(2 * np.pi * m * THETA) - a * np.sin(2 * np.pi * m * THETA)) * ( - np.sin(2 * np.pi * m * theta_c) + 2 * self.k * R / (m**2 * (2 * self.k * R + 1)) * ( - ((2 * np.pi * theta_c * m)**2 - 2) * np.sin(2 * np.pi * m * theta_c) + 4*np.pi*m*theta_c*np.cos(2 * np.pi * m * theta_c)))), axis=2) - - # Higher Fourier modes of change in power deficit wrt radius - dpdr += np.sum(((a * np.cos(2 * np.pi * m * THETA) + b * np.sin(2 * np.pi * m * THETA)) / (np.pi * m**3 * (2 * self.k * R + 1)**4) * ( - -4 * self.k * np.sin(2 * np.pi * m * theta_c) * (1 + m**2 + 2 * self.k * R * (m**2 - 2) + 2 * np.pi**2 * m**2 * (4 * self.k * R - 1) * theta_c**2) + - 2 * np.pi * m * np.cos(2 * np.pi * m * theta_c) * (4 * self.k * (1 - 4 * self.k * R) * theta_c + m**2 * (2 * self.k * R + 1) * ( - 1 + 2 * self.k * R + 8 * np.pi**2 * self.k * R * theta_c**2) * dtdr))), axis=2) - + dpdt = 0 + for m in np.arange(1,len(self.fs.b)): + # Higher Fourier modes of change in power deficit wrt angle + dpdt += (2 / (2 * self.k * R + 1)**2 * ( + self.fs.b[m] * np.cos(2 * np.pi * m * THETA) - self.fs.a[m] * np.sin(2 * np.pi * m * THETA)) * ( + np.sin(2 * np.pi * m * theta_c) + 2 * self.k * R / (m**2 * (2 * self.k * R + 1)) * ( + ((2 * np.pi * theta_c * m)**2 - 2) * np.sin(2 * np.pi * m * theta_c) + 4*np.pi*m*theta_c*np.cos(2 * np.pi * m * theta_c)))) + + # Higher Fourier modes of change in power deficit wrt radius + dpdr += ((self.fs.a[m] * np.cos(2 * np.pi * m * THETA) + self.fs.b[m] * np.sin(2 * np.pi * m * THETA)) / (np.pi * m**3 * (2 * self.k * R + 1)**4) * ( + -4 * self.k * np.sin(2 * np.pi * m * theta_c) * (1 + m**2 + 2 * self.k * R * (m**2 - 2) + 2 * np.pi**2 * m**2 * (4 * self.k * R - 1) * theta_c**2) + + 2 * np.pi * m * np.cos(2 * np.pi * m * theta_c) * (4 * self.k * (1 - 4 * self.k * R) * theta_c + m**2 * (2 * self.k * R + 1) * ( + 1 + 2 * self.k * R + 8 * np.pi**2 * self.k * R * theta_c**2) * dtdr))) + # Apply mask for points within rotor radius du = du * (1 - mask_area) + mask_val * mask_area np.fill_diagonal(du, 0.) @@ -162,8 +160,8 @@ def calculate_aep(self, gradient=False): # Complete gradient calculation if gradient==True: - dx = xx/np.sqrt(xx**2+yy**2)*dpdr + -yy/(2*np.pi*(xx**2+yy**2))*dpdt - dy = yy/np.sqrt(xx**2+yy**2)*dpdr + xx/(2*np.pi*(xx**2+yy**2))*dpdt + dx = xx/R*dpdr + -yy/(2*np.pi*R**2)*dpdt + dy = yy/R*dpdr + xx/(2*np.pi*R**2)*dpdt dx = np.nan_to_num(dx) dy = np.nan_to_num(dy) @@ -192,7 +190,7 @@ def calculate_aep(self, gradient=False): def _fourier_coefficients(self, num_terms=0): """ Compute the Fourier series expansion coefficients from the wind rose. - Modifies the Flowers interface in place to add a Fourier coefficients + Modifies the FlowersModel in place to add a Fourier coefficients dataframe: fs (pandas:dataframe): Fourier coefficients used to expand the wind rose: - 'a_free': real coefficients of freestream component @@ -206,7 +204,7 @@ def _fourier_coefficients(self, num_terms=0): """ # Resample wind rose for average wind speed per wind direction - wr = self.wind_rose.copy(deep=True) + wr = copy.deepcopy(self.wind_rose) wr = tl.resample_average_ws_by_wd(wr) # Transform wind direction to polar angle diff --git a/flowers/optimization/__init__.py b/flowers/optimization/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/model_interface.py b/flowers/optimization/model_interface.py similarity index 89% rename from model_interface.py rename to flowers/optimization/model_interface.py index ae85be8..793cbe8 100644 --- a/model_interface.py +++ b/flowers/optimization/model_interface.py @@ -1,12 +1,11 @@ import floris.tools as wfct import numpy as np -import pyoptsparse from scipy.interpolate import NearestNDInterpolator import time -import flowers_interface as flow -import optimization_interface as opt -import tools as tl +import flowers.flowers_model as flow +import flowers.optimization.optimization_interface as opt +import flowers.tools.tools as tl class AEPInterface(): @@ -25,7 +24,7 @@ class AEPInterface(): num_terms (int, optional): number of Fourier modes for the FLOWERS model k (float, optional): wake expansion rate in the FLOWERS model conventional_model (str, optional): underlying wake model: - - 'park' (default) + - 'jensen' (default) - 'gauss' turbine (str, optional): turbine type: - 'nrel_5MW' (default) @@ -45,10 +44,10 @@ def __init__(self, wind_rose, layout_x, layout_y, num_terms=0, k=0.05, conventio self.flowers_interface = flow.FlowersInterface(wind_rose, layout_x, layout_y, num_terms=num_terms, k=k, turbine=turbine) # Initialize FLORIS - if conventional_model is None or conventional_model == 'park': - self.floris_interface = wfct.floris_interface.FlorisInterface("./input/park.yaml") - elif conventional_model == 'gauss': - self.floris_interface = wfct.floris_interface.FlorisInterface("./input/gauss.yaml") + # if conventional_model is None or conventional_model == 'jensen': + # self.floris_interface = wfct.floris_interface.FlorisInterface("./input/jensen.yaml") + # elif conventional_model == 'gauss': + # self.floris_interface = wfct.floris_interface.FlorisInterface("./input/gauss.yaml") wd_array = np.array(wind_rose["wd"].unique(), dtype=float) ws_array = np.array(wind_rose["ws"].unique(), dtype=float) @@ -234,7 +233,7 @@ class WPLOInterface(): num_terms (int, optional): number of Fourier modes for the FLOWERS model k (float, optional): wake expansion rate in the FLOWERS model conventional_model (str, optional): underlying wake model: - - 'park' (default) + - 'jensen' (default) - 'gauss' turbine (str, optional): turbine type: - 'nrel_5MW' (default) @@ -250,8 +249,8 @@ def __init__(self, wind_rose, layout_x, layout_y, boundaries, num_terms=10, k=0. if conventional_model is None or conventional_model == 'gauss': self.floris_interface = wfct.floris_interface.FlorisInterface("./input/gauss.yaml") - elif conventional_model == 'park': - self.floris_interface = wfct.floris_interface.FlorisInterface("./input/park.yaml") + elif conventional_model == 'jensen': + self.floris_interface = wfct.floris_interface.FlorisInterface("./input/jensen.yaml") # Initialize FLOWERS interface self.flowers_interface = flow.FlowersInterface(wind_rose, layout_x, layout_y, num_terms=num_terms, k=k, turbine=turbine) @@ -265,7 +264,7 @@ def __init__(self, wind_rose, layout_x, layout_y, boundaries, num_terms=10, k=0. self.floris_interface.reinitialize(wind_directions=wr.wd,wind_speeds=wr.ws,layout_x=layout_x.flatten(),layout_y=layout_y.flatten(),time_series=True) # Initialize post-processing interface - self.post_processing = wfct.floris_interface.FlorisInterface("./input/gauss.yaml") + self.post_processing = wfct.floris_interface.FlorisInterface("./input/post.yaml") wind_rose = tl.resample_wind_speed(wind_rose, ws=np.arange(1.,26.,1.)) wd_array = np.array(wind_rose["wd"].unique(), dtype=float) ws_array = np.array(wind_rose["ws"].unique(), dtype=float) @@ -279,7 +278,7 @@ def __init__(self, wind_rose, layout_x, layout_y, boundaries, num_terms=10, k=0. self.post_processing.calculate_wake() self._aep_initial = np.sum(self.post_processing.get_farm_power() * self._freq_2D * 8760) - def run_optimization(self, optimizer, gradient="analytic", solver="SLSQP", timer=None, history='hist.hist', output='out.out'): + def run_optimization(self, optimizer, gradient="analytic", solver="SLSQP", scale=1e3, tol=1e-2, timer=None, history='hist.hist', output='out.out'): """ Run a Wind Plant Layout Optimization study with either the FLOWERS or Conventional optimizer. @@ -318,9 +317,34 @@ def run_optimization(self, optimizer, gradient="analytic", solver="SLSQP", timer # Instantiate optimizer class with user inputs if optimizer == "flowers": - prob = opt.FlowersOptimizer(self.flowers_interface, self._initial_x, self._initial_y, self._boundaries, grad=gradient, solver=solver, timer=timer, history_file=history, output_file=output) + prob = opt.FlowersOptimizer( + self.flowers_interface, + self._initial_x, + self._initial_y, + self._boundaries, + grad=gradient, + solver=solver, + scale=scale, + tol=tol, + timer=timer, + history_file=history, + output_file=output + ) elif optimizer == "conventional": - prob = opt.ConventionalOptimizer(self.floris_interface, self._freq_1D, self._initial_x, self._initial_y, self._boundaries, grad=gradient, solver=solver, timer=timer, history_file=history, output_file=output) + prob = opt.ConventionalOptimizer( + self.floris_interface, + self._freq_1D, + self._initial_x, + self._initial_y, + self._boundaries, + grad=gradient, + solver=solver, + scale=scale, + tol=tol, + timer=timer, + history_file=history, + output_file=output + ) # Solve optimization problem print("Solving layout optimization problem.") @@ -341,24 +365,17 @@ def run_optimization(self, optimizer, gradient="analytic", solver="SLSQP", timer self.solution["exit_code"] = sol.optInform['text'] self.solution["init_aep"] = self._aep_initial - # Get layout and objective history - hist = pyoptsparse.pyOpt_history.History(history,temp=True) - self.solution["hist_x"], self.solution["hist_y"] = prob.parse_hist_vars(hist) - self.solution["iter"] = len(self.solution["hist_x"]-1) - - # Post process AEP - hist_aep = np.zeros(len(self.solution["hist_x"])) - hist_aep[0] = self._aep_initial - for n in np.arange(1,self.solution["iter"]-1): - self.post_processing.reinitialize(layout_x=self.solution["hist_x"][n].flatten(),layout_y=self.solution["hist_y"][n].flatten()) - self.post_processing.calculate_wake() - hist_aep[n] = np.sum(self.post_processing.get_farm_power() * self._freq_2D * 8760) + # Get number of iterations + with open(output, 'r') as fp: + for line in fp: + if 'No. of major iterations' in line: + self.solution["iter"] = int(line.split()[4]) + break - self.post_processing.reinitialize(layout_x=self.solution["opt_x"].flatten(),layout_y=self.solution["opt_y"].flatten()) + # Compute optimal AEP + self.post_processing.reinitialize(layout_x=self.solution["opt_x"].flatten(),layout_y=self.solution["opt_y"].flatten(),time_series=False) self.post_processing.calculate_wake() self._aep_final = np.sum(self.post_processing.get_farm_power() * self._freq_2D * 8760) - hist_aep[-1] = self._aep_final self.solution["opt_aep"] = self._aep_final - self.solution["hist_aep"] = hist_aep return self.solution \ No newline at end of file diff --git a/optimization_interface.py b/flowers/optimization/optimization_interface.py similarity index 81% rename from optimization_interface.py rename to flowers/optimization/optimization_interface.py index 6982abb..e4cff21 100644 --- a/optimization_interface.py +++ b/flowers/optimization/optimization_interface.py @@ -9,7 +9,7 @@ class LayoutOptimizer(): """ - def _base_init_(self, layout_x, layout_y, boundaries, solver="SNOPT", timer=None, options=None, history_file='hist.hist', output_file='out.out'): + def _base_init_(self, layout_x, layout_y, boundaries, solver="SNOPT", tol=1e-2, timer=None, options=None, history_file='hist.hist', output_file='out.out'): # Save boundary information self._boundaries = np.array(boundaries).T @@ -27,13 +27,13 @@ def _base_init_(self, layout_x, layout_y, boundaries, solver="SNOPT", timer=None self._ymin = np.min(self._boundaries[1]) self._ymax = np.max(self._boundaries[1]) - self._x0 = self._norm(layout_x, self._xmin, self._xmax) - self._y0 = self._norm(layout_y, self._ymin, self._ymax) + self._x0 = layout_x + self._y0 = layout_y self._nturbs = len(layout_x) # Optimization initialization self.solver = solver - self.storeHistory = history_file + # self.storeHistory = history_file self.timeLimit = timer self.optProb = pyoptsparse.Optimization('layout', self._obj_func) @@ -46,13 +46,16 @@ def _base_init_(self, layout_x, layout_y, boundaries, solver="SNOPT", timer=None self.optOptions = options elif solver == "SNOPT": self.optOptions = { - "Print file": output_file, - "iSumm": 0, - "Major optimality tolerance": 1e-4, + # "Print file": output_file, + # "iSumm": 0, + "Summary file": output_file, + "iPrint": 0, + "Major optimality tolerance": tol, "Minor optimality tolerance": 1e-4, "Major feasibility tolerance": 1e-4, "Minor feasibility tolerance": 1e-4, "Scale option": 0, + # "Verify level": 3, } elif solver == "SLSQP": self.optOptions = { @@ -130,10 +133,14 @@ def _boundary_constraint(self, gradient=False): Cy[i] = (self._boundaries[0,idx] - self._boundaries[0,(idx+1)%self._nbounds]) / self._boundary_len[idx] # Distance is negative if inside boundary for optimization problem + # if gradient: + # return self._norm(C, self._xmin, self._xmax), Cx, Cy + # else: + # return self._norm(C, self._xmin, self._xmax) if gradient: - return self._norm(C, self._xmin, self._xmax), Cx, Cy + return C, Cx, Cy else: - return self._norm(C, self._xmin, self._xmax) + return C ########################################################################### # pyOptSparse wrapper functions @@ -141,25 +148,29 @@ def _boundary_constraint(self, gradient=False): def _optimize(self): if self.gradient: if self.timeLimit is not None: - self.sol = self.opt(self.optProb, storeHistory=self.storeHistory, sens=self._sens_func, timeLimit=self.timeLimit) + self.sol = self.opt(self.optProb, sens=self._sens_func)#, timeLimit=self.timeLimit) #storeHistory=self.storeHistory else: - self.sol = self.opt(self.optProb, storeHistory=self.storeHistory, sens=self._sens_func) + self.sol = self.opt(self.optProb, sens=self._sens_func) else: if self.timeLimit is not None: - self.sol = self.opt(self.optProb, storeHistory=self.storeHistory, timeLimit=self.timeLimit) + self.sol = self.opt(self.optProb, timeLimit=self.timeLimit) else: - self.sol = self.opt(self.optProb, storeHistory=self.storeHistory, sens='FD') + self.sol = self.opt(self.optProb, sens='FD') def parse_opt_vars(self, varDict): - self._x = self._unnorm(varDict["x"], self._xmin, self._xmax) - self._y = self._unnorm(varDict["y"], self._ymin, self._ymax) + # self._x = self._unnorm(varDict["x"], self._xmin, self._xmax) + # self._y = self._unnorm(varDict["y"], self._ymin, self._ymax) + self._x = varDict["x"] + self._y = varDict["y"] def parse_sol_vars(self, sol): - return np.array(self._unnorm(sol.getDVs()["x"], self._xmin, self._xmax)), np.array(self._unnorm(sol.getDVs()["y"], self._ymin, self._ymax)) + # return np.array(self._unnorm(sol.getDVs()["x"], self._xmin, self._xmax)), np.array(self._unnorm(sol.getDVs()["y"], self._ymin, self._ymax)) + return np.array(sol.getDVs()["x"]), np.array(sol.getDVs()["y"]) def parse_hist_vars(self, hist): val = hist.getValues(names=['x','y'],major=True) - return np.array(self._unnorm(val['x'], self._xmin, self._xmax)), np.array(self._unnorm(val['y'], self._ymin, self._ymax)) + # return np.array(self._unnorm(val['x'], self._xmin, self._xmax)), np.array(self._unnorm(val['y'], self._ymin, self._ymax)) + return np.array(val['x']), np.array(val['y']) def add_var_group(self, optProb): optProb.addVarGroup("x", self._nturbs, type="c", value=self._x0) @@ -193,15 +204,15 @@ class FlowersOptimizer(LayoutOptimizer): """ - def __init__(self, flowers_interface, layout_x, layout_y, boundaries, grad="analytical", solver="SNOPT", timer=None, history_file='hist.hist', output_file='out.out'): + def __init__(self, flowers_interface, layout_x, layout_y, boundaries, grad="analytical", solver="SNOPT", scale=1e3, tol=1e-2, timer=None, history_file='hist.hist', output_file='out.out'): self.model = flowers_interface if grad == "analytical": self.gradient = True elif grad == "numerical": self.gradient = False - self._aep_initial, self._grad_initial = self.model.calculate_aep(gradient=True) - self._grad_initial = 1e-4*self._aep_initial # TODO: adjust gradient scaling 1e-4 - self._base_init_(layout_x, layout_y, boundaries, solver=solver, timer=timer, history_file=history_file, output_file=output_file) + aep_initial = self.model.calculate_aep(gradient=False) + self._scale = aep_initial / scale + self._base_init_(layout_x, layout_y, boundaries, solver=solver, tol=tol, timer=timer, history_file=history_file, output_file=output_file) def _obj_func(self, varDict): # Parse the variable dictionary @@ -213,7 +224,7 @@ def _obj_func(self, varDict): # Compute the objective function funcs = {} funcs["obj"] = ( - -1 * self.model.calculate_aep() / self._aep_initial + -1 * self.model.calculate_aep() / self._scale ) # Compute constraints, if any are defined for the optimization @@ -231,7 +242,7 @@ def _sens_func(self, varDict, funcs): _, tmp = self.model.calculate_aep(gradient=True) funcsSens = {} - funcsSens["obj"] = {"x": -tmp[:,0]/self._grad_initial, "y": -tmp[:,1]/self._grad_initial} + funcsSens["obj"] = {"x": -tmp[:,0]/self._scale, "y": -tmp[:,1]/self._scale} _, tmpx, tmpy = self._boundary_constraint(gradient=True) funcsSens["con"] = {"x": np.diag(tmpx), "y": np.diag(tmpy)} @@ -259,12 +270,12 @@ class ConventionalOptimizer(LayoutOptimizer): """ - def __init__(self, floris_interface, freq_val, layout_x, layout_y, boundaries, grad="analytical", solver="SNOPT", timer=None, history_file='hist.hist', output_file='out.out'): + def __init__(self, floris_interface, freq_val, layout_x, layout_y, boundaries, grad="analytical", solver="SNOPT", scale=1e3, tol=1e-2, timer=None, history_file='hist.hist', output_file='out.out'): self.model = floris_interface self.gradient = False self._freq_1D = freq_val self.model.calculate_wake() - self._aep_initial = np.sum(self.model.get_farm_power() * self._freq_1D * 8760) + self._scale = np.sum(self.model.get_farm_power() * self._freq_1D * 8760) / scale self._base_init_(layout_x, layout_y, boundaries, solver=solver, timer=timer, history_file=history_file, output_file=output_file) def _obj_func(self, varDict): @@ -272,13 +283,13 @@ def _obj_func(self, varDict): self.parse_opt_vars(varDict) # Update turbine map with turbince locations - self.model.reinitialize(layout_x=self._x.flatten(), layout_y=self._y.flatten()) + self.model.reinitialize(layout_x=self._x.flatten(), layout_y=self._y.flatten(), time_series=True) # Compute the objective function self.model.calculate_wake() funcs = {} funcs["obj"] = ( - -1 * np.sum(self.model.get_farm_power() * self._freq_1D * 8760) / self._aep_initial + -1 * np.sum(self.model.get_farm_power() * self._freq_1D * 8760) / self._scale ) # Compute constraints, if any are defined for the optimization diff --git a/tools.py b/flowers/tools.py similarity index 79% rename from tools.py rename to flowers/tools.py index c2a81fa..b8cd234 100644 --- a/tools.py +++ b/flowers/tools.py @@ -1,12 +1,13 @@ -# FLOWERS -# Michael LoCascio +import copy import numpy as np import pandas as pd +import pickle from shapely.geometry import Polygon, Point -import floris.tools.wind_rose as rose + +# import floris.tools.wind_rose as rose ########################################################################### @@ -46,8 +47,8 @@ def random_layout(boundaries=[], n_turb=0, D=126.0, min_dist=2.0, idx=None): if idx != None: np.random.seed(idx) - xx = np.zeros(n_turb) - yy = np.zeros(n_turb) + xx = -1000*np.ones(n_turb) + yy = -1000*np.ones(n_turb) xmin = np.min([tup[0] for tup in boundaries]) xmax = np.max([tup[0] for tup in boundaries]) @@ -104,9 +105,10 @@ def discrete_layout(n_turb=0, D=126.0, min_dist=3.0, idx=None, spacing=False): yy = np.zeros(n_turb) # Indices of discrete grid - s = n_turb + 1 - x_idx = np.random.randint(0,s,n_turb) - y_idx = np.random.randint(0,s,n_turb) + sx = n_turb + 1 + sy = 6 + x_idx = np.random.randint(0,sx,n_turb) + y_idx = np.random.randint(0,sy,n_turb) pts = [(x_idx[i],y_idx[i]) for i in range(n_turb)] while len(np.unique(pts)) < len(pts): tmp = np.unique(pts) @@ -114,8 +116,8 @@ def discrete_layout(n_turb=0, D=126.0, min_dist=3.0, idx=None, spacing=False): for i in range(n_turb): if i not in tmp: new_set.append(i) - x_idx[new_set] = np.random.randint(0,s,len(new_set)) - y_idx[new_set] = np.random.randint(0,s,len(new_set)) + x_idx[new_set] = np.random.randint(0,sx,len(new_set)) + y_idx[new_set] = np.random.randint(0,sy,len(new_set)) pts = [(x_idx[i],y_idx[i]) for i in range(n_turb)] # Check that all combinations of x,y are unique @@ -133,102 +135,12 @@ def discrete_layout(n_turb=0, D=126.0, min_dist=3.0, idx=None, spacing=False): else: return xx, yy -def load_layout(name, boundaries=False): - """ - TODO: fill in description +def load_layout(idx, case, boundaries=True): + file = './layouts/' + case + str(idx) + '.p' + layout_x, layout_y, boundaries = pickle.load(open(file,'rb')) - """ - - if name == "iea": - layout_x = np.array([ - 2714.43, - 2416.08, - 1496.75, - 1860.65, - 2224.55, - 2588.45, - 1197.40, - 1619.09, - 2040.78, - 2462.46, - 898.05, - 1257.66, - 1617.27, - 1976.87, - 2336.48, - 598.70, - 1001.65, - 1404.60, - 1807.55, - 2210.50, - 299.35, - 750.21, - 1201.07, - 1651.93, - 2102.79, - 0.00, - 415.30, - 830.59, - 1245.88, - 1661.18, - 2076.47, - ]) - layout_y = np.array([ - 4042.95, - 3932.63, - 4020.72, - 3803.55, - 3586.38, - 3369.22, - 3618.82, - 3311.04, - 3003.27, - 2695.49, - 3216.92, - 2918.13, - 2619.34, - 2320.55, - 2021.76, - 2815.01, - 2448.27, - 2081.53, - 1714.78, - 1348.04, - 2413.11, - 1977.72, - 1542.33, - 1106.94, - 671.55, - 2011.21, - 1608.97, - 1206.72, - 804.48, - 402.24, - 0.00, - ]) - bound = [ - (2714.4, 4049.4), - (2132.7, 938.8), - (2092.8, 591.6), - (2078.9, 317.3), - (2076.1, 148.5), - (2076.6, 0.0), - (2076.5, 6.5), - (1208.6, 847.0), - (0.0, 2017.7), - (1496.7, 4027.2), - (1531.8, 4006.2), - (1931.2, 3818.5), - (2058.3, 3783.6), - (2192.8, 3792.9), - (2316.8, 3846.4), - (2416.0, 3939.1), - (2528.6, 4089.0), - (2550.9, 4126.3) - ] - if boundaries: - return layout_x, layout_y, bound + return layout_x, layout_y, boundaries else: return layout_x, layout_y @@ -236,34 +148,34 @@ def load_layout(name, boundaries=False): # Wind rose sampling ########################################################################### -def toolkit_wind_rose(lat, long): - """ - Sample a wind rose from the WIND Toolkit (copied from FLORIS) +# def toolkit_wind_rose(lat, long): +# """ +# Sample a wind rose from the WIND Toolkit (copied from FLORIS) - Args: - lat (float): latitude of wind farm site (in continental US) - long (float): longitude of wind farm site (in continental US) +# Args: +# lat (float): latitude of wind farm site (in continental US) +# long (float): longitude of wind farm site (in continental US) - Returns: - df (pandas.DataFrame): A dataframe for the wind rose in the FLORIS - format containing the following information: - - 'ws' (float): wind speeds [m/s] - - 'wd' (float): wind directions [deg] - - 'freq_val' (float): frequency for each wind speed and direction - - """ - - wind_rose = rose.WindRose() - wd_list = np.arange(0, 360, 1) - ws_list = np.arange(0, 26, 1) - df = wind_rose.import_from_wind_toolkit_hsds( - lat, - long, - ht=100, - wd=wd_list, - ws=ws_list, - ) - return df +# Returns: +# df (pandas.DataFrame): A dataframe for the wind rose in the FLORIS +# format containing the following information: +# - 'ws' (float): wind speeds [m/s] +# - 'wd' (float): wind directions [deg] +# - 'freq_val' (float): frequency for each wind speed and direction + +# """ + +# wind_rose = rose.WindRose() +# wd_list = np.arange(0, 360, 1) +# ws_list = np.arange(0, 26, 1) +# df = wind_rose.import_from_wind_toolkit_hsds( +# lat, +# long, +# ht=100, +# wd=wd_list, +# ws=ws_list, +# ) +# return df def load_wind_rose(idx): """ @@ -423,13 +335,13 @@ def resample_wind_speed(df, ws=np.arange(0, 26, 1.0)): return df -def resample_average_ws_by_wd(df): +def resample_average_ws_by_wd(wind_rose): """ Calculate the mean wind speed for each wind direction bin and resample the wind rose. (Copied from FLORIS) Args: - df (pandas.DataFrame): Wind rose DataFrame containing the following + wind_rose (WindRose): Wind rose object containing the following columns: - 'wd': Wind direction bin center values (deg). - 'ws': Wind speed bin center values (m/s). @@ -442,10 +354,11 @@ def resample_average_ws_by_wd(df): - 'ws': Resampled average wind speed values (m/s). - 'freq_val': The resampled frequency of occurance of the wind conditions in the other columns. - + """ - # Make a copy of incoming dataframe - df = df.copy(deep=True) + # Make a copy of incoming FLORIS WindRose object + df = copy.deepcopy(wind_rose) + # wind_rose = copy.deepcopy(wind_rose) ws_avg = [] @@ -457,6 +370,20 @@ def resample_average_ws_by_wd(df): / df.loc[df["wd"] == val]["freq_val"].sum() ) + # freq_avg = np.sum(wind_rose.freq_table, axis=1) + # ws_avg = np.sum( + # wind_rose.freq_table * wind_rose.wind_speeds.reshape(1, -1), axis=1 + # ) / freq_avg + + # print(wind_rose.wind_directions) + # print(wind_rose.freq_table) + # print(ws_avg) + # lkj + # df = pd.DataFrame() + # df["wd"] = wind_rose.wind_directions + # df["ws"] = ws_avg + # df["freq_val"] = freq_avg + # Regroup df = df.groupby("wd").sum() @@ -482,10 +409,10 @@ def ct_lookup(u, turbine_type, ct=None): Args: u (float): normalized inflow wind speed - + Returns: ct (float): thrust coefficient - + """ if ct != None: @@ -503,7 +430,7 @@ def ct_lookup(u, turbine_type, ct=None): 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5, 12.0, 12.5, 13.0, 13.5, 14.0, 14.5, 15.0, 15.5, 16.0, 16.5, 17.0, 17.5, 18.0, 18.5, 19.0, 19.5, 20.0, 20.5, 21.0, 21.5, 22.0, 22.5, 23.0, 23.5, 24.0, 24.5, 25.0, 25.01, 25.02, 50.0]) - + return np.interp(u, u_table, ct_table) def cp_lookup(u, turbine_type, cp=None): @@ -512,10 +439,10 @@ def cp_lookup(u, turbine_type, cp=None): Args: u (float): normalized inflow wind speed - + Returns: cp (float): power coefficient - + """ if cp != None: cp_table = np.array([0.0, 0.0, cp, cp, 0.0, 0.0]) diff --git a/visualization.py b/flowers/visualization.py similarity index 92% rename from visualization.py rename to flowers/visualization.py index ea4fdf6..162ab21 100644 --- a/visualization.py +++ b/flowers/visualization.py @@ -7,7 +7,7 @@ import matplotlib.pyplot as plt import numpy as np -import tools as tl +import flowers.tools.tools as tl ########################################################################### @@ -76,6 +76,8 @@ def plot_wind_rose( width=0.9 * np.radians(wd_step), color=color_array(ws_idx), edgecolor="k", + linewidth=0.5, + label=ws_labels[ws_idx], ) ) # break @@ -91,8 +93,10 @@ def plot_wind_rose( ax.set_theta_direction(-1) ax.set_theta_offset(np.pi / 2.0) ax.set_theta_zero_location("N") + ax.tick_params(axis='x', which='major', pad=-1) # 1 for AEP paper ax.set_xticklabels(["N", "NE", "E", "SE", "S", "SW", "W", "NW"]) ax.set_yticklabels([]) + ax.set_axisbelow(True) return ax @@ -219,31 +223,45 @@ def plot_optimal_layout( y0 = y_init y1 = y_final verts = boundaries - r = D/2 + r = D xlab = 'x [m]' ylab = 'y [m]' # Plot turbine locations - ax.scatter(x0, y0, s=0.01, color=color_initial) + # ax.scatter(x0, y0, s=0.01, color=color_initial) ax.scatter(x1, y1, s=0.01, color=color_final) - for x, y in zip(x0, y0): - ax.add_patch(plt.Circle((x, y), r, color=color_initial)) + # for x, y in zip(x0, y0): + # ax.add_patch(plt.Circle((x, y), r, color=color_initial)) for x, y in zip(x1, y1): ax.add_patch(plt.Circle((x, y), r, color=color_final)) ax.set(xlabel=xlab, ylabel=ylab, aspect='equal') - ax.legend(['Initial','Final'],markerscale=50) + ax.legend(['Final'],markerscale=50) + # ax.legend(['Initial','Final'],markerscale=50) # leg.legendHandles[0]._legmarker.set_markersize(6) # leg.legendHandles[1]._legmarker.set_markersize(6) ax.grid() - # Plot plant boundary - for i in range(len(verts)): - if i == len(verts) - 1: - ax.plot([verts[i][0], verts[0][0]], [verts[i][1], verts[0][1]], "black") - else: - ax.plot( - [verts[i][0], verts[i + 1][0]], [verts[i][1], verts[i + 1][1]], "black" - ) + if verts.ndim == 2: + # Plot plant boundary for single boundary + for i in range(len(verts)): + if i == len(verts) - 1: + ax.plot([verts[i][0], verts[0][0]], [verts[i][1], verts[0][1]], "black") + else: + ax.plot( + [verts[i][0], verts[i + 1][0]], [verts[i][1], verts[i + 1][1]], "black" + ) + + else: + # Plot plant boundary for multiple boundaries + for ii in range(len(verts)): + sub_verts = verts[ii] + for i in range(len(sub_verts)): + if i == len(sub_verts) - 1: + ax.plot([sub_verts[i][0], sub_verts[0][0]], [sub_verts[i][1], sub_verts[0][1]], "black") + else: + ax.plot( + [sub_verts[i][0], sub_verts[i + 1][0]], [sub_verts[i][1], sub_verts[i + 1][1]], "black" + ) return ax diff --git a/generate_wind_roses.py b/generate_wind_roses.py deleted file mode 100644 index a951e3b..0000000 --- a/generate_wind_roses.py +++ /dev/null @@ -1,33 +0,0 @@ -import matplotlib.pyplot as plt -import os -import numpy as np -import pandas as pd - -import visualization as vis - -wd = np.arange(0., 360., 1.0) -ws = 10. -a = 1 -c = 2 -# freq = a*(np.exp(-(wd-270)**2/(2*c)**2)) -freq = np.ones(len(wd)) -freq[270] = 1e7/20 -freq /= np.sum(freq) - -df = pd.DataFrame() -df["wd"] = wd -df["ws"] = ws -df["freq_val"] = freq - -vis.plot_wind_rose(df) - -plt.show(block=False) - -var = input("Would you like to save this wind rose? [y/n]: ") - -if var == 'y': - idx = input("Provide an index for this wind rose: ") - file_name = './wind_roses/wr' + idx + '.p' - df.to_pickle(file_name) - -plt.close() \ No newline at end of file diff --git a/gradient_debug.py b/gradient_debug.py deleted file mode 100644 index f26da97..0000000 --- a/gradient_debug.py +++ /dev/null @@ -1,134 +0,0 @@ -import flowers_interface as flow -import numpy as np -import tools as tl -import visualization as viz -import matplotlib.pyplot as plt -import old_files.model as set -from scipy.interpolate import NearestNDInterpolator -import floris.tools as wfct - -import warnings -warnings.filterwarnings("ignore") - -wr = tl.load_wind_rose(6) # 14 -# wr["wd"] = np.remainder(450 - wr.wd, 360) -# wr.sort_values("wd", inplace=True) - -# layout_x, layout_y = tl.discrete_layout(6, idx=12) - -# layout_x = 126. * np.array([0.,7.,14.]) -# layout_y = 126. * np.array([0.,0.,0.]) -# layout_x = 126. * np.array([0.,0.,7.,7.]) -# layout_y = 126. * np.array([0.,7.,0.,7.]) -layout_x = 126. * np.array([0.,0.,0.,7.,7.,7.,14.,14.,14.]) -layout_y = 126. * np.array([0.,7.,14.,0.,7.,14.,0.,7.,14.]) - - -fi = flow.FlowersInterface(wr,layout_x,layout_y) -aep, grad = fi.calculate_aep(gradient=True) - -# grad = fi.calculate_aep_gradient() - - -# fli = wfct.floris_interface.FlorisInterface('./input/park.yaml') -# fli.reinitialize( -# layout_x=layout_x.flatten(), -# layout_y=layout_y.flatten(), -# wind_shear=0) -# # Initialize wind direction-speed frequency array for AEP -# wd_array = np.array(wr["wd"].unique(), dtype=float) -# ws_array = np.array(wr["ws"].unique(), dtype=float) -# wd_grid, ws_grid = np.meshgrid(wd_array, ws_array, indexing="ij") -# freq_interp = NearestNDInterpolator(wr[["wd", "ws"]], wr["freq_val"]) -# freq = freq_interp(wd_grid, ws_grid) -# freq_floris = freq / np.sum(freq) - -# fli.reinitialize( -# wind_directions=wd_array, -# wind_speeds=ws_array) -# fli.calculate_wake() -# aep_floris = np.sum(fli.get_farm_power() * freq_floris * 8760) -# print(aep_floris) -# dasd - - -# fi.reinitialize(num_terms=10) -# aep, grad_new = fi.calculate_aep(gradient=True) -# print(grad) -# print(grad_new) -# print((grad-grad_new)/grad*100) -# dasd - -fd = np.zeros((len(layout_x),2)) -h = 126.*1e-8 - -for i in range(len(layout_x)): - layout_xx = np.copy(layout_x) - layout_xx[i] += h - fi.reinitialize(layout_x=layout_xx) - fd1 = fi.calculate_aep() - layout_xx[i] -= 2*h - fi.reinitialize(layout_x=layout_xx) - fd2 = fi.calculate_aep() - - fd[i,0] = (fd1 - fd2)/(2*h) - # fd[i,0] = (fd1 - aep)/h - - -for i in range(len(layout_x)): - layout_yy = np.copy(layout_y) - layout_yy[i] += h - fi.reinitialize(layout_y=layout_yy) - fd1 = fi.calculate_aep() - layout_yy[i] -= 2*h - fi.reinitialize(layout_y=layout_yy) - fd2 = fi.calculate_aep() - - fd[i,1] = (fd1 - fd2)/(2*h) - # fd[i,1] = (fd1 - aep)/h - -# fd_new = np.zeros((len(layout_x),2)) -# fi.fourier_coefficients(num_terms=5) -# aep_new = fi.calculate_aep() - -# for i in range(len(layout_x)): -# layout_xx = np.copy(layout_x) -# layout_xx[i] += h -# tmp = flow.Flowers(wr,layout_xx,layout_y) -# tmp.fourier_coefficients(num_terms=5) -# fd1 = tmp.calculate_aep() - -# layout_xx[i] -= 2*h -# tmp = flow.Flowers(wr,layout_xx,layout_y) -# tmp.fourier_coefficients(num_terms=5) -# fd2 = tmp.calculate_aep() - -# fd_new[i,0] = (fd1 - fd2)/(2*h) - -# for i in range(len(layout_x)): -# layout_yy = np.copy(layout_y) -# layout_yy[i] += h -# tmp = flow.Flowers(wr,layout_x,layout_yy) -# tmp.fourier_coefficients(num_terms=5) -# fd1 = tmp.calculate_aep() - -# layout_yy[i] -= 2*h -# tmp = flow.Flowers(wr,layout_x,layout_yy) -# tmp.fourier_coefficients(num_terms=5) -# fd2 = tmp.calculate_aep() - -# fd_new[i,1] = (fd1 - fd2)/(2*h) - -# print(fd/aep*100) -# print(fd_new/aep_new*100) -# print((fd_new-fd)/fd*100) - -print(grad) -print(fd) -# print(grad/aep*100) -print(np.abs((grad-fd)/fd)*100) -print(np.sign(grad/fd)) - -viz.plot_wind_rose(wr) -viz.plot_layout(layout_x,layout_y) -plt.show() \ No newline at end of file diff --git a/input/gauss.yaml b/input/gauss.yaml deleted file mode 100644 index 5035bb8..0000000 --- a/input/gauss.yaml +++ /dev/null @@ -1,89 +0,0 @@ - -name: Gauss-Jimenez -description: Three turbines using Jensen / Jimenez models -floris_version: v3.0.0 - -logging: - console: - enable: true - level: WARNING - file: - enable: false - level: WARNING - -solver: - type: turbine_grid - turbine_grid_points: 1 - -farm: - layout_x: - - 0.0 - - 630.0 - - 1260.0 - layout_y: - - 0.0 - - 0.0 - - 0.0 - turbine_type: - - nrel_5MW - -flow_field: - air_density: 1.225 - reference_wind_height: -1 # -1 is code for use the hub height - turbulence_intensity: 0.07 - wind_directions: - - 270.0 - wind_shear: 0.0 - wind_speeds: - - 8.0 - wind_veer: 0.0 - -wake: - model_strings: - combination_model: fls - deflection_model: gauss - turbulence_model: crespo_hernandez - velocity_model: gauss - - enable_secondary_steering: false - enable_yaw_added_recovery: false - enable_transverse_velocities: false - - wake_deflection_parameters: - gauss: - ad: 0.0 - alpha: 0.58 - bd: 0.0 - beta: 0.077 - dm: 1.0 - ka: 0.38 - kb: 0.004 - jimenez: - ad: 0.0 - bd: 0.0 - kd: 0.05 - - wake_velocity_parameters: - cc: - a_s: 0.179367259 - b_s: 0.0118889215 - c_s1: 0.0563691592 - c_s2: 0.13290157 - a_f: 3.11 - b_f: -0.68 - c_f: 2.41 - alpha_mod: 1.0 - gauss: - alpha: 0.58 - beta: 0.077 - ka: 0.38 - kb: 0.004 - jensen: - we: 0.05 - - wake_turbulence_parameters: - crespo_hernandez: - initial: 0.1 - constant: 0.5 - ai: 0.8 - downstream: -0.32 diff --git a/input/park.yaml b/input/park.yaml deleted file mode 100644 index 62b61a1..0000000 --- a/input/park.yaml +++ /dev/null @@ -1,89 +0,0 @@ - -name: Park -description: Three turbines using Jensen / Jimenez models -floris_version: v3.0.0 - -logging: - console: - enable: true - level: WARNING - file: - enable: false - level: WARNING - -solver: - type: turbine_grid - turbine_grid_points: 1 - -farm: - layout_x: - - 0.0 - - 630.0 - - 1260.0 - layout_y: - - 0.0 - - 0.0 - - 0.0 - turbine_type: - - nrel_5MW - -flow_field: - air_density: 1.225 - reference_wind_height: -1 # -1 is code for use the hub height - turbulence_intensity: 0.06 - wind_directions: - - 270.0 - wind_shear: 0.0 - wind_speeds: - - 8.0 - wind_veer: 0.0 - -wake: - model_strings: - combination_model: fls - deflection_model: jimenez - turbulence_model: crespo_hernandez - velocity_model: jensen - - enable_secondary_steering: false - enable_yaw_added_recovery: false - enable_transverse_velocities: false - - wake_deflection_parameters: - gauss: - ad: 0.0 - alpha: 0.58 - bd: 0.0 - beta: 0.077 - dm: 1.0 - ka: 0.38 - kb: 0.004 - jimenez: - ad: 0.0 - bd: 0.0 - kd: 0.05 - - wake_velocity_parameters: - cc: - a_s: 0.179367259 - b_s: 0.0118889215 - c_s1: 0.0563691592 - c_s2: 0.13290157 - a_f: 3.11 - b_f: -0.68 - c_f: 2.41 - alpha_mod: 1.0 - gauss: - alpha: 0.58 - beta: 0.077 - ka: 0.38 - kb: 0.004 - jensen: - we: 0.05 - - wake_turbulence_parameters: - crespo_hernandez: - initial: 0.1 - constant: 0.5 - ai: 0.8 - downstream: -0.32 diff --git a/opt_debug.py b/opt_debug.py deleted file mode 100644 index d1a3ec7..0000000 --- a/opt_debug.py +++ /dev/null @@ -1,60 +0,0 @@ -import model_interface as inter -import tools as tl -import numpy as np -import visualization as vis -import matplotlib.pyplot as plt -import matplotlib.animation as anim - -wr = tl.load_wind_rose(8) -layout_x = 126. * np.array([0.,0.,-3.,7.,7.,7.,14.,14.,14.]) -layout_y = 126. * np.array([0.,7.,14.,0.,7.,14.,-2.,7.,14.]) -boundaries = [(0., 0.),(12*126, 0.),(12*126, 6*126.),(6*126, 12*126.),(0, 12*126.)] - -# wr = tl.load_wind_rose(1) -# layout_x = 126. * np.array([0.,0.,1.,7.,7.,7.,14.,14.,14.,18.,18.,18.]) -# layout_y = 126. * np.array([0.,7.,9.,0.,7.,10.,0.,7.,14.,-1.,5.,10.]) -# boundaries = [(0., 0.),(20*126, 0.),(20*126, 20*126.),(15*126, 20*126.),(15*126, 10*126.),(0, 10*126.)] - -# layout_x = 126. * np.array([-2.,-2.,-2.,5.,5.,5.,12.,12.,12.]) -# layout_y = 126. * np.array([-2.,5.,12.,-2.,5.,12.,-2.,5.,12.]) -# boundaries = [(0., 0.),(10*126, 0.),(10*126, 10*126.),(0, 10*126.)] - -opt = inter.WPLOInterface(wr, layout_x, layout_y, boundaries) -solution = opt.run_optimization(optimizer="conventional", solver="SNOPT", gradient="numerical", timer=60) - -print(solution["init_aep"]) -print(solution["opt_aep"]) -print(solution["total_time"]) -print(solution["iter"]) -print(solution["obj_calls"]) -print("Exit code: " + str(solution["exit_code"])) - -vis.plot_optimal_layout(np.array(boundaries), solution["opt_x"],solution["opt_y"],solution["init_x"],solution["init_y"]) -# plt.savefig('./figures/opt_example.png', dpi=500) -vis.plot_wind_rose(wr) - -plt.figure() -plt.plot(range(solution["iter"]),solution["hist_aep"]/1e9) -plt.xlabel('Iteration') -plt.ylabel('AEP [GWh]') -# plt.savefig('./figures/opt_conv_example.png', dpi=500) - -# fig, ax = plt.subplots(1,1) -# ax.set(aspect='equal', xlim=[-5,25], ylim=[-5,25], xlabel='x/D', ylabel='y/D') -# # line, = ax.plot(solution["hist_x"][0]/126.,solution["hist_y"][0]/126.,"o",color='tab:red',markersize=7) -# patches = [] -# for n in range(len(solution["hist_x"][0])): -# patches.append(ax.add_patch(plt.Circle((solution["hist_x"][0,n], solution["hist_y"][0,n]), 1/2, color='tab:red'))) - -# def animate(i): -# for n in range(len(solution["hist_x"][0])): -# patches[n].center = solution["hist_x"][i,n]/126., solution["hist_y"][i,n]/126. -# # line.set_data(solution["hist_x"][i]/126., solution["hist_y"][i]/126.) -# ax.set(title=str(i)) - -# mov = anim.FuncAnimation(fig, animate, frames=solution["iter"], repeat=True) - -# mov.save("./figures/opt_example.gif", dpi=500) - -plt.show() - diff --git a/opt_figures.py b/opt_figures.py deleted file mode 100644 index ee4aa71..0000000 --- a/opt_figures.py +++ /dev/null @@ -1,292 +0,0 @@ -import numpy as np -import tools as tl -import matplotlib.cm as cm -import matplotlib.pyplot as plt -import matplotlib.colors as co -import pickle -import visualization as vis -from scipy.optimize import curve_fit -import matplotlib.animation as animation -import matplotlib.collections as coll -from mpl_toolkits.axes_grid1 import make_axes_locatable -from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, mark_inset - -save = True -dpi = 500 -figs = "parameter" - -########################################################################### -# FLOWERS MULTISTART -########################################################################### - -if figs == "multistart": - # Load files - N = 25 # idx = 8 outlier - solutions = [] - for n in range(N): - file_name = 'solutions/opt_multi' + str(n) + '.p' - solution, wr, boundaries = pickle.load(open(file_name,'rb')) - solutions.append(solution) - - boundaries /= 126. - boundaries = np.append(boundaries,np.reshape(boundaries[:,0],(2,1)),axis=1) - - vis.plot_wind_rose(wr) - - fig, ax = plt.subplots(1,1) - # aep_cmap = cm.get_cmap('winter') - # nn = np.linspace(0.,1.,N,endpoint=True) - for n in range(N): - solution = solutions[n] - ax.plot(range(solution['iter']),solution['hist_aep']/1e9,color='tab:blue',alpha=0.7) - ax.set(xlabel='Iteration', ylabel='AEP [GWh]') - fig.tight_layout() - if save: - plt.savefig('./figures/opt_multistart_history.png', dpi=dpi) - - fig, ax = plt.subplots(1,1) - for n in range(N): - solution = solutions[n] - layout = [] - for i in range(len(solution['opt_x'])): - layout.append(plt.Circle((solution['opt_x'][i]/126., solution['opt_y'][i]/126.), 1/2)) - - layouts = coll.PatchCollection(layout, color='tab:blue', alpha=0.5) - ax.add_collection(layouts) - - ax.plot(boundaries[0],boundaries[1],color='k',linewidth=2,zorder=1) - ax.set(xlabel='x/D', ylabel='y/D', aspect='equal', title='Optimal Layouts') - ax.grid(True) - fig.tight_layout() - if save: - plt.savefig('./figures/opt_multistart_opt_layouts.png', dpi=dpi) - - fig, ax = plt.subplots(1,1) - for n in range(N): - solution = solutions[n] - layout = [] - for i in range(len(solution['init_x'])): - layout.append(plt.Circle((solution['init_x'][i]/126., solution['init_y'][i]/126.), 1/2)) - - layouts = coll.PatchCollection(layout, color='tab:orange', alpha=0.5) - ax.add_collection(layouts) - - ax.plot(boundaries[0],boundaries[1],color='k',linewidth=2,zorder=1) - ax.set(xlabel='x/D', ylabel='y/D', aspect='equal', title='Initial Layouts') - ax.grid(True) - fig.tight_layout() - if save: - plt.savefig('./figures/opt_multistart_init_layouts.png', dpi=dpi) - - # 16 is optimal case - fig, ax = plt.subplots(1,1) - nn = 16 - solution = solutions[nn] - layout_init = [] - layout_final = [] - for i in range(len(solution['opt_x'])): - layout_init.append(plt.Circle((solution['init_x'][i]/126., solution['init_y'][i]/126.), 1/2)) - layout_final.append(plt.Circle((solution['opt_x'][i]/126., solution['opt_y'][i]/126.), 1/2)) - - tmp0 = plt.Circle(([],[]),1/2,color='tab:orange',label='Initial Layout') - tmp1 = plt.Circle(([],[]),1/2,color='tab:blue',label='Optimal Layout') - - layout0 = coll.PatchCollection(layout_init, color='tab:orange') - layout1 = coll.PatchCollection(layout_final, color='tab:blue') - ax.add_collection(layout0) - ax.add_collection(layout1) - ax.plot(boundaries[0],boundaries[1],color='k',linewidth=2,zorder=1) - ax.set(xlabel='x/D', ylabel='y/D', aspect='equal',title='Case ' + str(nn)) - ax.grid(True) - ax.legend(handles=[tmp0,tmp1],loc='upper right') - fig.tight_layout() - if save: - plt.savefig('./figures/opt_multistart_success.png', dpi=dpi) - - fig, ax = plt.subplots(1,1) - for n in range(N): - solution = solutions[n] - ax.scatter(solution['total_time'],solution['opt_aep']/1e9,20,color='tab:blue') - ax.set(xlabel='Solver Time [s]',ylabel='Optimal AEP [GWh]') - # print(ax.get_ylim()) - if save: - plt.savefig('./figures/opt_multistart_aepvstime.png', dpi=dpi) - - fig, ax = plt.subplots(1,1) - for n in range(N): - solution = solutions[n] - ax.scatter(solution['iter'],solution['opt_aep']/1e9,20,color='tab:blue') - ax.set(xlabel='Iterations',ylabel='Optimal AEP [GWh]') - if save: - plt.savefig('./figures/opt_multistart_aepvsiter.png', dpi=dpi) - - fig, ax = plt.subplots(1,1) - for n in range(N): - solution = solutions[n] - ax.scatter(solution['init_aep']/1e9,solution['opt_aep']/1e9,20,color='tab:blue') - ax.set(xlabel='Initial AEP [GWh]',ylabel='Optimal AEP [GWh]') - xlim = ax.get_xlim() - ylim = ax.get_ylim() - xmin = np.min([xlim[0],ylim[0]]) - xmax = np.max([xlim[1],ylim[1]]) - ax.set(xlim=[xmin,xmax], ylim=[xmin,xmax]) - if save: - plt.savefig('./figures/opt_multistart_aepvsaep.png', dpi=dpi) - - fig, ax = plt.subplots(1,1) - for n in range(N): - solution = solutions[n] - ax.scatter(solution['iter'],solution['total_time'],20,color='tab:blue') - ax.set(xlabel='Iterations',ylabel='Solver Time [s]') - if save: - plt.savefig('./figures/opt_multistart_timevsiter.png', dpi=dpi) - - fig, ax = plt.subplots(1,1,figsize=(13,3)) - exit_codes = {} - for n in range(N): - exit_code = solutions[n]['exit_code'] - if exit_code in exit_codes: - exit_codes[exit_code] += 1 - else: - exit_codes[exit_code] = 1 - - ax.bar(range(len(exit_codes)), list(exit_codes.values()), align='center',color='tab:blue',width=0.2) - ax.set(xlabel='SNOPT Exit Codes', ylabel='Count',xticks=range(len(exit_codes)), xticklabels=list(exit_codes.keys())) - fig.tight_layout() - if save: - plt.savefig('./figures/opt_multistart_codes.png', dpi=dpi) - - optimal_aep = np.zeros(N) - solver_time = np.zeros(N) - for n in range(N): - optimal_aep[n] = solutions[n]['opt_aep'] - solver_time[n] = solutions[n]['total_time'] - - print("Median Optimal AEP: {:.2f} GWh".format(np.median(optimal_aep)/1e9)) - print("Median Solver Time: {:.2f} s".format(np.median(solver_time))) - - plt.show() - -########################################################################### -# FLOWERS PARAMETER SENSITIVITY -########################################################################### - -if figs == "parameter": - # Load files - k = np.array([0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.10]) - solutions = [] - - for idx in range(len(k)): - file_name = 'solutions/opt_k' + str(k[idx])[-2:] + '.p' - solution, wr, boundaries = pickle.load(open(file_name,'rb')) - solutions.append(solution) - - boundaries /= 126. - boundaries = np.append(boundaries,np.reshape(boundaries[:,0],(2,1)),axis=1) - - - fig, ax = plt.subplots(1,1) - for n in range(len(k)): - solution = solutions[n] - ax.plot(range(solution['iter']),solution['hist_aep']/1e9,color='tab:blue',alpha=0.7) - ax.set(xlabel='Iteration', ylabel='AEP [GWh]') - fig.tight_layout() - if save: - plt.savefig('./figures/opt_parameter_history.png', dpi=dpi) - - fig, ax = plt.subplots(1,1) - for n in range(len(k)): - solution = solutions[n] - layout = [] - # if n == 1: - # continue - for i in range(len(solution['opt_x'])): - layout.append(plt.Circle((solution['opt_x'][i]/126., solution['opt_y'][i]/126.), 1/2)) - - layouts = coll.PatchCollection(layout, color='tab:blue', alpha=0.5) - ax.add_collection(layouts) - - ax.plot(boundaries[0],boundaries[1],color='k',linewidth=2,zorder=1) - ax.set(xlabel='x/D', ylabel='y/D', aspect='equal', title='Optimal Layouts') - ax.grid(True) - fig.tight_layout() - if save: - plt.savefig('./figures/opt_parameter_layouts.png', dpi=dpi) - - fig, ax = plt.subplots(2,5,figsize=(14,7)) - for nn in range(len(k)): - axx = int(nn >= 5) - axy = int(nn % 5) - solution = solutions[nn] - layout_init = [] - layout_final = [] - for i in range(len(solution['opt_x'])): - layout_init.append(plt.Circle((solution['init_x'][i]/126., solution['init_y'][i]/126.), 1/2)) - layout_final.append(plt.Circle((solution['opt_x'][i]/126., solution['opt_y'][i]/126.), 1/2)) - - tmp0 = plt.Circle(([],[]),1/2,color='tab:orange',label='Initial Layout') - tmp1 = plt.Circle(([],[]),1/2,color='tab:blue',label='Optimal Layout') - - layout0 = coll.PatchCollection(layout_init, color='tab:orange') - layout1 = coll.PatchCollection(layout_final, color='tab:blue') - ax[axx,axy].add_collection(layout0) - ax[axx,axy].add_collection(layout1) - ax[axx,axy].plot(boundaries[0],boundaries[1],color='k',linewidth=2,zorder=1) - ax[axx,axy].set(xlabel='x/D', ylabel='y/D', aspect='equal',title='$k$ = ' + str(k[nn])) - ax[axx,axy].grid(True) - # ax[axx,axy].legend(handles=[tmp0,tmp1],loc='upper right') - fig.tight_layout() - if save: - plt.savefig('./figures/opt_parameter_layoutsidx.png', dpi=dpi) - - fig, ax = plt.subplots(1,1) - for n in range(len(k)): - solution = solutions[n] - ax.scatter(k[n],solution['opt_aep']/1e9,20,color='tab:blue') - # xlim = ax.get_xlim() - # ax.hlines(solution['init_aep']/1e9,xlim[0],xlim[1],color='tab:blue',linestyles='--') - ax.set(xlabel='$k$',ylabel='Optimal AEP [GWh]', ylim=[388.60,392.09]) - if save: - plt.savefig('./figures/opt_parameter_aep.png', dpi=dpi) - - fig, ax = plt.subplots(1,1) - for n in range(len(k)): - solution = solutions[n] - ax.scatter(k[n],solution['total_time'],20,color='tab:blue') - ax.set(xlabel='$k$',ylabel='Solver Time [s]') - if save: - plt.savefig('./figures/opt_parameter_time.png', dpi=dpi) - - fig, ax = plt.subplots(1,1) - for n in range(len(k)): - solution = solutions[n] - ax.scatter(k[n],solution['iter'],20,color='tab:blue') - ax.set(xlabel='$k$',ylabel='Iterations') - if save: - plt.savefig('./figures/opt_parameter_iter.png', dpi=dpi) - - fig, ax = plt.subplots(1,1,figsize=(13,3)) - exit_codes = {} - for n in range(len(k)): - exit_code = solutions[n]['exit_code'] - if exit_code in exit_codes: - exit_codes[exit_code] += 1 - else: - exit_codes[exit_code] = 1 - - ax.bar(range(len(exit_codes)), list(exit_codes.values()), align='center',color='tab:blue',width=0.2) - ax.set(xlabel='SNOPT Exit Codes', ylabel='Count',xticks=range(len(exit_codes)), xticklabels=list(exit_codes.keys())) - fig.tight_layout() - if save: - plt.savefig('./figures/opt_parameter_codes.png', dpi=dpi) - - # optimal_aep = np.zeros(len(k)) - # solver_time = np.zeros(len(k)) - # for n in range(2): - # optimal_aep[n] = solutions[n]['opt_aep'] - # solver_time[n] = solutions[n]['total_time'] - - # print("Median Optimal AEP: {:.2f} GWh".format(np.median(optimal_aep)/1e9)) - # print("Median Solver Time: {:.2f} s".format(np.median(solver_time))) - - plt.show() \ No newline at end of file diff --git a/opt_multistart.py b/opt_multistart.py deleted file mode 100644 index f50e651..0000000 --- a/opt_multistart.py +++ /dev/null @@ -1,29 +0,0 @@ -import model_interface as inter -import numpy as np -import tools as tl -import pickle - -idx = 0 -model = "conventional" -gradients = "numerical" - -wr = tl.load_wind_rose(1) -boundaries = [(3*126., 0.),(12*126, 0.),(15*126, 10*126.),(10*126, 15*126.),(0, 15*126.)] -file_base = 'solutions/opt_floris' - -save_file = file_base + str(idx) + '.p' -history_file = file_base + str(idx) + '.hist' -output_file = file_base + str(idx) + '.out' - -layout_x, layout_y = tl.random_layout(boundaries, n_turb=20, idx=idx) -opt = inter.WPLOInterface(wr, layout_x, layout_y, boundaries) - -solution = opt.run_optimization(optimizer=model, solver="SNOPT", gradient=gradients, timer=600, history=history_file, output=output_file) - -boundaries = np.array(boundaries).T - -print(solution['init_aep']) -print(solution['opt_aep']) -print(solution['total_time']) -print(solution['iter']) -# pickle.dump((solution, wr, boundaries), open(save_file,'wb')) \ No newline at end of file diff --git a/opt_parameter.py b/opt_parameter.py deleted file mode 100644 index 806a391..0000000 --- a/opt_parameter.py +++ /dev/null @@ -1,26 +0,0 @@ -import model_interface as inter -import numpy as np -import tools as tl -import pickle - -idx = 16 -model = "flowers" -gradients = "analytical" -# k = [0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10] -k = 0.10 - -wr = tl.load_wind_rose(1) -boundaries = [(3*126., 0.),(12*126, 0.),(15*126, 10*126.),(10*126, 15*126.),(0, 15*126.)] -file_base = 'solutions/opt_k' - -save_file = file_base + str(k)[-2:] + '.p' -history_file = file_base + str(k)[-2:] + '.hist' -output_file = file_base + str(k)[-2:] + '.out' - -layout_x, layout_y = tl.random_layout(boundaries, n_turb=20, idx=idx) -opt = inter.WPLOInterface(wr, layout_x, layout_y, boundaries, k=k) - -solution = opt.run_optimization(optimizer=model, solver="SNOPT", gradient=gradients, timer=300, history=history_file, output=output_file) - -boundaries = np.array(boundaries).T -pickle.dump((solution, wr, boundaries), open(save_file,'wb')) \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..8e7b84b --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,58 @@ + +# Project configuration for FLOWERS +# TOML spec: https://toml.io/en/v1.0.0 + +[build-system] +requires = [ + "setuptools", +] +build-backend = "setuptools.build_meta" + +# https://setuptools.pypa.io/en/latest/userguide/package_discovery.html +[tool.setuptools] +packages = ["flowers"] + +[project] +name = "FLOWERS" +version = "1.0.0-alpha0" +authors = [ + {name = "Christopher Bay", email = "christopher.bay@nrel.gov"}, +] +description = "A package for rapid energy simulation of wind farms" +readme= "README.md" +requires-python = ">=3.10, <3.13" +classifiers = [ + "Development Status :: 3 - Alpha", + "License :: OSI Approved :: Apache Software License", + "Operating System :: OS Independent", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + "Programming Language :: Python :: 3.12", + # "Programming Language :: Python :: 3.13", +] +dependencies = [ + "numpy", + "matplotlib", + "scipy", + "pandas", + "shapely", +] +[project.optional-dependencies] +dev = [ + "black", + "pytest", + "pytest-cov", +] +docs = [ + "pyxdsm", + "jupyter-book", + "sphinx-book-theme", + "sphinx-autodoc-typehints", +] +[project.urls] +# Homepage = "https://example.com" +Documentation = "https://wisdem.github.io/flowers" +Repository = "https://github.com/wisdem/flowers.git" +Issues = "https://github.com/wisdem/flowers/issues" +# Changelog = "https://github.com/wisdem/flowers/blob/main/CHANGELOG.md" \ No newline at end of file diff --git a/save_wind_roses.py b/save_wind_roses.py deleted file mode 100644 index 56eb2d8..0000000 --- a/save_wind_roses.py +++ /dev/null @@ -1,30 +0,0 @@ -# FLOWERS - -# Michael LoCascio - -import matplotlib.pyplot as plt - -import tools as tl -import visualization as vis - -""" -This file downloads a wind rose from the WIND Toolkit and gives the -user the option to save it to a file '.wind_rose/wr#.p', where # is a -user-provided index. - -""" - -wind_rose = tl.toolkit_wind_rose(lat = 41.09, long = -72.2) - -vis.plot_wind_rose(wind_rose) - -plt.show(block=False) - -var = input("Would you like to save this wind rose? [y/n]: ") - -if var == 'y': - idx = input("Provide an index for this wind rose: ") - file_name = './wind_roses/wr' + idx + '.p' - wind_rose.to_pickle(file_name) - -plt.close() \ No newline at end of file diff --git a/show_wind_roses.py b/show_wind_roses.py deleted file mode 100644 index 5e05d88..0000000 --- a/show_wind_roses.py +++ /dev/null @@ -1,23 +0,0 @@ -# FLOWERS - -# Michael LoCascio - -import matplotlib.pyplot as plt -import os -import pandas as pd - -import visualization as vis - -""" -This file visualizes all of the wind roses saved in the -./wind_rose/ directory. - -""" - -for file in os.listdir('./wind_roses/'): - if file[-1] == 'p': - df = pd.read_pickle('./wind_roses/' + file) - vis.plot_wind_rose(df) - ax = plt.gca() - ax.set(title=file) -plt.show() \ No newline at end of file diff --git a/solutions/.DS_Store b/solutions/.DS_Store deleted file mode 100644 index 5008ddf..0000000 Binary files a/solutions/.DS_Store and /dev/null differ diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/regression/test_flowers_aep.py b/tests/regression/test_flowers_aep.py new file mode 100644 index 0000000..3d531f6 --- /dev/null +++ b/tests/regression/test_flowers_aep.py @@ -0,0 +1,34 @@ + +from pathlib import Path + +import numpy as np +import pandas as pd +from pytest import approx + +import flowers +from flowers import FlowersModel + + +def test_regression_aep(): + # Create generic layout + D = 126. + layout_x = D * np.array([0.,0.,0.,7.,7.,7.,14.,14.,14.,21.,21.,21.,28.,28.,28.]) + layout_y = D * np.array([0.,7.,14.,0.,7.,14.,0.,7.,14.,0.,7.,14.,0.,7.,14.]) + + # Load in wind data + wind_rose_file = Path( + Path(flowers.__file__).parent, + "..", + "examples", + "inputs", + "HKW_wind_rose.csv", + ) + df = pd.read_csv(wind_rose_file) + + # Setup FLOWERS model + flowers_model = FlowersModel(df, layout_x, layout_y) + + # Calculate the AEP + aep = flowers_model.calculate_aep() + + assert aep == approx(350958996034.8524, 1e-3) \ No newline at end of file diff --git a/wind_roses/.DS_Store b/wind_roses/.DS_Store deleted file mode 100644 index 5008ddf..0000000 Binary files a/wind_roses/.DS_Store and /dev/null differ diff --git a/wind_roses/wr0.p b/wind_roses/wr0.p deleted file mode 100644 index b2ec851..0000000 Binary files a/wind_roses/wr0.p and /dev/null differ diff --git a/wind_roses/wr00.p b/wind_roses/wr00.p deleted file mode 100644 index 6fad596..0000000 Binary files a/wind_roses/wr00.p and /dev/null differ diff --git a/wind_roses/wr1.p b/wind_roses/wr1.p deleted file mode 100644 index 93ea276..0000000 Binary files a/wind_roses/wr1.p and /dev/null differ diff --git a/wind_roses/wr10.p b/wind_roses/wr10.p deleted file mode 100644 index f8a5c85..0000000 Binary files a/wind_roses/wr10.p and /dev/null differ diff --git a/wind_roses/wr11.p b/wind_roses/wr11.p deleted file mode 100644 index c63f9fa..0000000 Binary files a/wind_roses/wr11.p and /dev/null differ diff --git a/wind_roses/wr12.p b/wind_roses/wr12.p deleted file mode 100644 index de49985..0000000 Binary files a/wind_roses/wr12.p and /dev/null differ diff --git a/wind_roses/wr13.p b/wind_roses/wr13.p deleted file mode 100644 index 6fad596..0000000 Binary files a/wind_roses/wr13.p and /dev/null differ diff --git a/wind_roses/wr14.p b/wind_roses/wr14.p deleted file mode 100644 index 2c81f04..0000000 Binary files a/wind_roses/wr14.p and /dev/null differ diff --git a/wind_roses/wr15.p b/wind_roses/wr15.p deleted file mode 100644 index 79a7e1a..0000000 Binary files a/wind_roses/wr15.p and /dev/null differ diff --git a/wind_roses/wr2.p b/wind_roses/wr2.p deleted file mode 100644 index 0a7e4f5..0000000 Binary files a/wind_roses/wr2.p and /dev/null differ diff --git a/wind_roses/wr3.p b/wind_roses/wr3.p deleted file mode 100644 index 88f258f..0000000 Binary files a/wind_roses/wr3.p and /dev/null differ diff --git a/wind_roses/wr4.p b/wind_roses/wr4.p deleted file mode 100644 index 29259ab..0000000 Binary files a/wind_roses/wr4.p and /dev/null differ diff --git a/wind_roses/wr5.p b/wind_roses/wr5.p deleted file mode 100644 index ccea68b..0000000 Binary files a/wind_roses/wr5.p and /dev/null differ diff --git a/wind_roses/wr6.p b/wind_roses/wr6.p deleted file mode 100644 index 873bdcf..0000000 Binary files a/wind_roses/wr6.p and /dev/null differ diff --git a/wind_roses/wr7.p b/wind_roses/wr7.p deleted file mode 100644 index 8a1c897..0000000 Binary files a/wind_roses/wr7.p and /dev/null differ diff --git a/wind_roses/wr8.p b/wind_roses/wr8.p deleted file mode 100644 index 50ca3e2..0000000 Binary files a/wind_roses/wr8.p and /dev/null differ diff --git a/wind_roses/wr9.p b/wind_roses/wr9.p deleted file mode 100644 index 83648a9..0000000 Binary files a/wind_roses/wr9.p and /dev/null differ