Iterate over param values in biosnicar v2

biosnicar_iterator.py

#!/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=",")