Created
June 2, 2023 08:41
-
-
Save jmcook1186/b8b494a83a5736c79e55593ab7d82e2e to your computer and use it in GitHub Desktop.
Iterate over param values in biosnicar v2
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/env python3 | |
# -*- coding: utf-8 -*- | |
# save this file to the top level biosnicar_py directory, then run from the command line | |
# as follows: | |
# python src/biosnicar/biosnicar_iterator.py | |
import numpy as np | |
from pathlib import Path | |
from adding_doubling_solver import adding_doubling_solver | |
from column_OPs import get_layer_OPs, mix_in_impurities | |
from display import display_out_data, plot_albedo | |
from setup_snicar import setup_snicar | |
from toon_rt_solver import toon_solver | |
from validate_inputs import validate_inputs | |
BIOSNICAR_SRC_PATH = Path(__file__).resolve().parent | |
# define input file | |
input_file = BIOSNICAR_SRC_PATH.joinpath("inputs.yaml").as_posix() | |
# run setup_snicar() to establish base-case values for ALL necessary params | |
( | |
ice, | |
illumination, | |
rt_config, | |
model_config, | |
plot_config, | |
impurities, | |
) = setup_snicar(input_file) | |
# now define the range of values you actually want to iterate over | |
lyrList = [0, 1] | |
densList = [400, 500, 600, 700, 800] | |
reffList = [200, 400, 600, 800, 1000] | |
zenList = [30, 40, 50, 60] | |
bcList = [500, 1000, 2000] | |
dzList = [ | |
[0.08, 0.1], | |
[0.10, 0.15], | |
[0.2, 0.5], | |
[0.3, 0.5], | |
[1, 10], | |
] | |
# set up an output array | |
ncols = ( | |
len(lyrList) | |
* len(densList) | |
* len(reffList) | |
* len(zenList) | |
* len(bcList) | |
* len(dzList) | |
) | |
specOut = np.zeros(shape=(ncols, 481)) | |
# iterate over all your values | |
counter = 0 | |
for layer_type in lyrList: | |
for density in densList: | |
for reff in reffList: | |
for zen in zenList: | |
for bc in bcList: | |
for dz in dzList: | |
ice.dz = dz | |
ice.nbr_lyr = 2 | |
ice.layer_type = [layer_type] * len(ice.dz) | |
ice.rho = [density] * len(ice.dz) | |
ice.rds = [reff] * len(ice.dz) | |
illumination.solzen = zen | |
# remember to recalculate irradiance values when any of the dependency values change | |
# i.e. irradiance is derived from solzen, so call the recalc func after updating solzen | |
illumination.calculate_irradiance() | |
impurities[0].conc = [ | |
bc, | |
bc, | |
] # bc in all layers | |
# remember to recalculate RI after changing any ice optical values | |
ice.calculate_refractive_index(input_file) | |
ssa_snw, g_snw, mac_snw = get_layer_OPs(ice, model_config) | |
tau, ssa, g, L_snw = mix_in_impurities( | |
ssa_snw, | |
g_snw, | |
mac_snw, | |
ice, | |
impurities, | |
model_config, | |
) | |
# now call the solver of your choice (here AD solver) | |
outputs = adding_doubling_solver( | |
tau, ssa, g, L_snw, ice, illumination, model_config | |
) | |
# spectral albedo appended to output array | |
specOut[counter, 0:480] = outputs.albedo | |
specOut[counter, 480] = outputs.BBA | |
counter += 1 | |
# save to file | |
np.savetxt("/home/joe/Desktop/py_benchmark_data.csv", specOut, delimiter=",") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment