text:
import DataFrames.DataFrame
import DataArrays.skipmissing

include("clusters.jl")

sim_param = setup_sim_param_model()
add_param_fixed(sim_param,"num_targets_sim_pass_one",150061) #9)   # For "observed" data, use a realistic number of targets (after any cuts you want to perform)
add_param_fixed(sim_param,"max_incl_sys",80.0) #degrees; 0 (deg) for isotropic system inclinations; set closer to 90 (deg) for more transiting systems

const max_incl_sys = get_real(sim_param,"max_incl_sys")
cos_factor = cos(max_incl_sys*pi/180) #factor to divide the number of targets in simulation by to get the actual number of targets needed (with an isotropic distribution of system inclinations) to produce as many transiting systems for a single observer

##### To load the DR25 catalog and compute arrays of multiplicities, periods, period ratios, transit durations, transit depths, period-normalized transit duration ratios (xi), and transit depth ratios:

Q1Q17_DR25 = CSV.read("q1_q17_dr25_koi.tab_selectcols_new.csv", header=19, nullable=true)
Q1Q17_DR25_stellar = CSV.read("q1_q17_dr25_stellar_koi.tab_selectcols.csv", header=30, nullable=true)
Q1Q17_DR25_stellar_all = CSV.read("q1_q17_dr25_stellar_koi.tab_all.csv", header=1, nullable=true)

N_Kepler_targets = sum((Q1Q17_DR25_stellar_all[:teff] .> 4000.) .& (Q1Q17_DR25_stellar_all[:teff] .< 7000.) .& (Q1Q17_DR25_stellar_all[:logg] .> 4.))
println("Total number of Kepler targets satisfying our cuts: ", N_Kepler_targets)

table_confirmed = Q1Q17_DR25[(Q1Q17_DR25[:koi_disposition] .== "CONFIRMED") .| (Q1Q17_DR25[:koi_disposition] .== "CANDIDATE"), :] #Table containing only the confirmed and candidate objects
table_stellar = Q1Q17_DR25_stellar

table_confirmed = table_confirmed[(table_confirmed[:koi_period] .> get_real(sim_param,"min_period")) .& (table_confirmed[:koi_period] .< get_real(sim_param,"max_period")), :] #to make additional cuts in period P to be comparable to our simulated sample
table_confirmed = table_confirmed[(table_confirmed[:koi_prad] .> get_real(sim_param,"min_radius")/ExoplanetsSysSim.earth_radius) .& (table_confirmed[:koi_prad] .< get_real(sim_param,"max_radius")/ExoplanetsSysSim.earth_radius) .& (.~ismissing.(table_confirmed[:koi_prad])), :] #to make additional cuts in planetary radii to be comparable to our simulated sample

#To make cuts based on stellar properties of T_eff and logg:
teff_confirmed = zeros(Int64, size(table_confirmed,1)) #list to be filled with the T_eff (K) for the objects
logg_confirmed = zeros(Float64, size(table_confirmed,1)) #list to be filled with the logg(cgs) for the objects
cdpp5_confirmed = zeros(Float64, size(table_confirmed,1)) #list to be filled with the RMS CDPP 5h values for the objects
cdpp6_confirmed = zeros(Float64, size(table_confirmed,1)) #list to be filled with the RMS CDPP 6h values for the objects
for (i,KepID) in enumerate(table_confirmed[:kepid])
    teff_confirmed[i] = table_stellar[:teff][table_stellar[:kepid] .== KepID, :][1]
    logg_confirmed[i] = table_stellar[:logg][table_stellar[:kepid] .== KepID, :][1]
    cdpp5_confirmed[i] = table_stellar[:rrmscdpp05p0][table_stellar[:kepid] .== KepID, :][1]
    cdpp6_confirmed[i] = table_stellar[:rrmscdpp06p0][table_stellar[:kepid] .== KepID, :][1]
end

cdpp_cut = 250.
println("Fraction of CONFIRMED and CANDIDATE planets after CDPP cut: ", size(table_confirmed[(teff_confirmed .> 4000.) .& (teff_confirmed .< 7000.) .& (logg_confirmed .> 4.) .& (cdpp5_confirmed .< cdpp_cut), :], 1), "/", size(table_confirmed[(teff_confirmed .> 4000.) .& (teff_confirmed .< 7000.) .& (logg_confirmed .> 4.), :], 1))

table_confirmed = table_confirmed[(teff_confirmed .> 4000.) .& (teff_confirmed .< 7000.) .& (logg_confirmed .> 4.) .& (cdpp5_confirmed .< cdpp_cut), :]

KOI_systems = [x[1:6] for x in table_confirmed[:kepoi_name]]
checked_bools = zeros(size(table_confirmed,1)) #0's denote KOI that were not checked yet; 1's denote already checked KOI

M_confirmed = Int64[] #list to be filled with the planet multiplicities of the systems
R_confirmed = Float64[] #list to be filled with period ratios of adjacent planet pairs
xi_confirmed = Float64[] #list to be filled with the period-normalized transit duration ratios of adjacent planet pairs
D_ratio_confirmed = Float64[] #list to be filled with the transit depth ratios of adjacent planet pairs
P_confirmed = table_confirmed[:koi_period] #array of the periods (days)
P_confirmed = collect(skipmissing(P_confirmed))
t_D_confirmed = table_confirmed[:koi_duration] #array of the transit durations (hrs)
t_D_confirmed = collect(skipmissing(t_D_confirmed))
D_confirmed = table_confirmed[:koi_depth]/(1e6) #array of the transit depths (fraction)
D_confirmed = D_confirmed[ismissing.(D_confirmed) .== false] #to get rid of NA values

for i in 1:length(KOI_systems)
    if checked_bools[i] == 0 #if the KOI has not been checked (included while looking at another planet in the same system)
        system_i = (1:length(KOI_systems))[KOI_systems .== KOI_systems[i]]
        checked_bools[system_i] = 1

        #To get the periods and transit durations in this system:
        system_P = table_confirmed[:koi_period][system_i] #periods of all the planets in this system
        system_t_D = table_confirmed[:koi_duration][system_i] #transit durations of all the planets in this system
        system_D = table_confirmed[:koi_depth][system_i] #transit depths (in ppm) of all the planets in this system
        system_sort_i = sortperm(system_P) #indices that would sort the periods of the planets in this system
        system_P = system_P[system_sort_i] #periods of all the planets in this system, sorted
        system_t_D = system_t_D[system_sort_i] #transit durations of all the planets in this system, sorted by period
        system_D = system_D[system_sort_i] #transit depths of all the planets in this system, sorted by period

        #To count the total number of planets in this system:
        push!(M_confirmed, length(system_P))

        #To compute the period ratios, period-normalized transit duration ratios, and transit depth ratios in this system:
        system_R = system_P[2:end]./system_P[1:end-1] #period ratios of all the adjacent planet pairs in this system
        system_D_ratio = system_D[2:end]./system_D[1:end-1] #transit depth ratios of all the adjacent planet pairs in this system
        system_xi = (system_t_D[1:end-1]./system_t_D[2:end]).*(system_P[2:end]./system_P[1:end-1]).^(1//3) #period-normalized transit duration ratios of all the adjacent planet pairs in this system

        append!(R_confirmed, system_R)
        append!(D_ratio_confirmed, system_D_ratio)
        append!(xi_confirmed, system_xi)
    end
end

##### To run a set of simulations on a grid of two parameters, saving KS statistics:

#=
#To run our model on a 2D grid of 'power_law_P' and 'power_law_r':

power_law_P_range = linspace(-0.5, 0.5, 11)
power_law_r_range = linspace(-3., -1.2, 11)

model_fit_stats = zeros(6,length(power_law_P_range),length(power_law_r_range)) #3D array to be filled with all the statistics of the model fits: number of observed planets, and KS distances for multiplicities, periods, period ratios, transit depths, and period-normalized transit duration ratios
model_fit_stats_fields = ["Number of observed planets", "KS multiplicities", "KS periods", "KS period ratios", "KS transit depths", "KS period-normalized transit duration ratios"] #names of the fields associated with 'model_fit_stats' array

sim_param_default = setup_sim_param_model() #To set up all the default parameters

for (i,power_law_P) in enumerate(power_law_P_range)
    for (j,power_law_r) in enumerate(power_law_r_range)
        tic()
        sim_param_here = deepcopy(sim_param_default) #To set up all the default parameters
        add_param_fixed(sim_param_here,"power_law_P",power_law_P) #To replace the 'power_law_P' parameter
        add_param_fixed(sim_param_here,"power_law_r",power_law_r) #To replace the 'power_law_r' parameter

        cat_phys = generate_kepler_physical_catalog(sim_param_here)
        cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys,sim_param_here)
        cat_obs = observe_kepler_targets_single_obs(cat_phys_cut,sim_param_here)
        summary_stat = calc_summary_stats_model(cat_obs,sim_param_here)

        M_cat_obs = ones(Int64,0) #array to be filled with the number of transiting planets in each simulated system
        for k in 1:get_int(sim_param_here,"max_tranets_in_sys")
            append!(M_cat_obs, k*ones(Int64, summary_stat.stat["num n-tranet systems"][k]))
        end

        model_fit_stats[1,i,j] = summary_stat.stat["num_tranets"]
        model_fit_stats[2,i,j] = ksstats_ints(M_cat_obs, M_confirmed)[5]
        model_fit_stats[3,i,j] = ksstats(summary_stat.stat["P list"], P_confirmed)[5]
        model_fit_stats[4,i,j] = ksstats(summary_stat.stat["period_ratio_list"], R_confirmed)[5]
        model_fit_stats[5,i,j] = ksstats(summary_stat.stat["depth list"], D_confirmed)[5]
        model_fit_stats[6,i,j] = ksstats(summary_stat.stat["duration_ratio_list"], xi_confirmed)[5]

        println("# power_law_P, power_law_r: ", get_real(sim_param_here,"power_law_P"), ", ",  get_real(sim_param_here,"power_law_r"))
        toc()
    end
end
=#

#=
f = open("Model_stats_P_r_grid.out", "w")
println(f,"# Model statistics on a grid of (power_law_P, power_law_r) = (", power_law_P_range, ", ",  power_law_r_range, ")")
println(f,"# num_targets_sim_pass_one: ", get_int(sim_param_default,"num_targets_sim_pass_one"))
println(f,"# log_rate_clusters: ", get_real(sim_param_default,"log_rate_clusters"))
println(f,"# log_rate_planets_per_cluster: ", get_real(sim_param_default,"log_rate_planets_per_cluster"))
println(f,"# sigma_incl: ", get_real(sim_param_default,"sigma_incl"))
println(f,"# sigma_incl_near_mmr: ", get_real(sim_param_default,"sigma_incl_near_mmr"))
println(f,"# sigma_hk: ", get_real(sim_param_default,"sigma_hk"))
println(f,"# num_mutual_hill_radii: ", get_real(sim_param_default,"num_mutual_hill_radii"))
println(f,"# mr_power_index: ", get_real(sim_param_default,"mr_power_index"))
println(f,"# mr_max_mass: ", get_real(sim_param_default,"mr_max_mass"))
println(f,"# sigma_logperiod_per_pl_in_cluster: ", get_real(sim_param_default,"sigma_logperiod_per_pl_in_cluster"))
for i in 1:size(model_fit_stats)[1]
    println(f,"# ", model_fit_stats_fields[i], ":")
    println(f,model_fit_stats[i,:,:]')
end
close(f)
=#

#=
#To run our model on a 2D grid of 'num_mutual_hill_radii' and 'sigma_logperiod_per_pl_in_cluster':

num_mutual_hill_radii_range = linspace(5., 20., 11)
sigma_logperiod_per_pl_in_cluster_range = linspace(0.05, 0.5, 11)

model_fit_stats = zeros(8,length(num_mutual_hill_radii_range),length(sigma_logperiod_per_pl_in_cluster_range)) #3D array to be filled with all the statistics of the model fits: number of observed planets, and KS distances for multiplicities, periods, period ratios, transit durations, period-normalized transit duration ratios, transit depths, and transit depth ratios
model_fit_stats_fields = ["Number of observed planets", "KS multiplicities", "KS periods", "KS period ratios", "KS transit durations", "KS period-normalized transit duration ratios", "KS transit depths", "KS transit depth ratios"] #names of the fields associated with 'model_fit_stats' array

for (i,num_mutual_hill_radii) in enumerate(num_mutual_hill_radii_range)
    for (j,sigma_logperiod_per_pl_in_cluster) in enumerate(sigma_logperiod_per_pl_in_cluster_range)
        tic()
        sim_param_here = deepcopy(sim_param) #To set up all the default parameters
        add_param_fixed(sim_param_here,"num_mutual_hill_radii",num_mutual_hill_radii) #To replace the 'num_mutual_hill_radii' parameter
        add_param_fixed(sim_param_here,"sigma_logperiod_per_pl_in_cluster",sigma_logperiod_per_pl_in_cluster) #To replace the 'sigma_logperiod_per_pl_in_cluster' parameter

        cat_phys = generate_kepler_physical_catalog(sim_param_here)
        cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys,sim_param_here)
        cat_obs = observe_kepler_targets_single_obs(cat_phys_cut,sim_param_here)
        summary_stat = calc_summary_stats_model(cat_obs,sim_param_here)

        M_cat_obs = ones(Int64,0) #array to be filled with the number of transiting planets in each simulated system
        for k in 1:get_int(sim_param_here,"max_tranets_in_sys")
            append!(M_cat_obs, k*ones(Int64, summary_stat.stat["num n-tranet systems"][k]))
        end

        model_fit_stats[1,i,j] = summary_stat.stat["num_tranets"]
        model_fit_stats[2,i,j] = ksstats_ints(M_cat_obs, M_confirmed)[5]
        model_fit_stats[3,i,j] = ksstats(summary_stat.stat["P list"], P_confirmed)[5]
        model_fit_stats[4,i,j] = ksstats(summary_stat.stat["period_ratio_list"], R_confirmed)[5]
        model_fit_stats[5,i,j] = ksstats(summary_stat.stat["duration list"].*24, t_D_confirmed)[5] #transit durations in simulations are in days, while in the Kepler catalog are in hours
        model_fit_stats[6,i,j] = ksstats(summary_stat.stat["duration_ratio_list"], xi_confirmed)[5]
        model_fit_stats[7,i,j] = ksstats(summary_stat.stat["depth list"], D_confirmed)[5]
        model_fit_stats[8,i,j] = ksstats(summary_stat.stat["radius_ratio_list"].^2, D_ratio_confirmed)[5] #simulations save radius ratios while we computed transit duration ratios from the Kepler catalog

        println("# num_mutual_hill_radii, sigma_logperiod_per_pl_in_cluster: ", get_real(sim_param_here,"num_mutual_hill_radii"), ", ",  get_real(sim_param_here,"sigma_logperiod_per_pl_in_cluster"))
        toc()
    end
end
=#

#=
f = open("Model_stats_Nhill_sigmaP_grid.out", "w")
println(f,"# Model statistics on a grid of (num_mutual_hill_radii, sigma_logperiod_per_pl_in_cluster) = (", num_mutual_hill_radii_range, ", ",  sigma_logperiod_per_pl_in_cluster_range, ")")
println(f, "# num_targets_sim_pass_one: ", get_int(sim_param,"num_targets_sim_pass_one"))
println(f, "# max_incl_sys: ", get_real(sim_param,"max_incl_sys"))
println(f, "# log_rate_clusters: ", get_real(sim_param,"log_rate_clusters"))
println(f, "# max_clusters_in_sys: ", get_int(sim_param,"max_clusters_in_sys"))
println(f, "# log_rate_planets_per_cluster: ", get_real(sim_param,"log_rate_planets_per_cluster"))
println(f, "# max_planets_in_clusters: ", get_int(sim_param,"max_planets_in_cluster"))
println(f, "# power_law_P: ", get_real(sim_param,"power_law_P"))
println(f, "# min_period: ", get_real(sim_param,"min_period"))
println(f, "# max_period: ", get_real(sim_param,"max_period"))
println(f, "# power_law_r1: ", get_real(sim_param,"power_law_r1"))
println(f, "# power_law_r2: ", get_real(sim_param,"power_law_r2"))
println(f, "# min_radius (R_earth): ", get_real(sim_param,"min_radius")/ExoplanetsSysSim.earth_radius)
println(f, "# max_radius (R_earth): ", get_real(sim_param,"max_radius")/ExoplanetsSysSim.earth_radius)
println(f, "# break_radius (R_earth): ", get_real(sim_param,"break_radius")/ExoplanetsSysSim.earth_radius)
println(f, "# sigma_incl: ", get_real(sim_param,"sigma_incl"))
println(f, "# sigma_incl_near_mmr: ", get_real(sim_param,"sigma_incl_near_mmr"))
println(f, "# sigma_hk: ", get_real(sim_param,"sigma_hk"))
#####println(f, "# num_mutual_hill_radii: ", get_real(sim_param,"num_mutual_hill_radii"))
println(f, "# mr_power_index: ", get_real(sim_param,"mr_power_index"))
println(f, "# mr_max_mass (M_earth): ", get_real(sim_param,"mr_max_mass")/ExoplanetsSysSim.earth_mass)
println(f, "# sigma_log_radius_in_cluster: ", get_real(sim_param,"sigma_log_radius_in_cluster"))
#####println(f, "# sigma_logperiod_per_pl_in_cluster: ", get_real(sim_param,"sigma_logperiod_per_pl_in_cluster"))
for i in 1:size(model_fit_stats)[1]
    println(f,"# ", model_fit_stats_fields[i], ":")
    println(f,model_fit_stats[i,:,:]')
end
close(f)
=#

#=
#To run our model on a 2D grid of 'sigma_incl' and 'sigma_hk':

sigma_incl_range = linspace(0., 5., 11)
sigma_hk_range = linspace(0., 0.2, 11)

model_fit_stats = zeros(8,length(sigma_incl_range),length(sigma_hk_range)) #3D array to be filled with all the statistics of the model fits: number of observed planets, and KS distances for multiplicities, periods, period ratios, transit durations, period-normalized transit duration ratios, transit depths, and transit depth ratios
model_fit_stats_fields = ["Number of observed planets", "KS multiplicities", "KS periods", "KS period ratios", "KS transit durations", "KS period-normalized transit duration ratios", "KS transit depths", "KS transit depth ratios"] #names of the fields associated with 'model_fit_stats' array

for (i,sigma_incl) in enumerate(sigma_incl_range)
    for (j,sigma_hk) in enumerate(sigma_hk_range)
        tic()
        sim_param_here = deepcopy(sim_param) #To set up all the default parameters
        add_param_fixed(sim_param_here,"sigma_incl",sigma_incl) #To replace the 'sigma_incl' parameter
        add_param_fixed(sim_param_here,"sigma_hk",sigma_hk) #To replace the 'sigma_hk' parameter

        cat_phys = generate_kepler_physical_catalog(sim_param_here)
        cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys,sim_param_here)
        cat_obs = observe_kepler_targets_single_obs(cat_phys_cut,sim_param_here)
        summary_stat = calc_summary_stats_model(cat_obs,sim_param_here)

        M_cat_obs = ones(Int64,0) #array to be filled with the number of transiting planets in each simulated system
        for k in 1:get_int(sim_param_here,"max_tranets_in_sys")
            append!(M_cat_obs, k*ones(Int64, summary_stat.stat["num n-tranet systems"][k]))
        end

        model_fit_stats[1,i,j] = summary_stat.stat["num_tranets"]
        model_fit_stats[2,i,j] = ksstats_ints(M_cat_obs, M_confirmed)[5]
        model_fit_stats[3,i,j] = ksstats(summary_stat.stat["P list"], P_confirmed)[5]
        model_fit_stats[4,i,j] = ksstats(summary_stat.stat["period_ratio_list"], R_confirmed)[5]
        model_fit_stats[5,i,j] = ksstats(summary_stat.stat["duration list"].*24, t_D_confirmed)[5] #transit durations in simulations are in days, while in the Kepler catalog are in hours
        model_fit_stats[6,i,j] = ksstats(summary_stat.stat["duration_ratio_list"], xi_confirmed)[5]
        model_fit_stats[7,i,j] = ksstats(summary_stat.stat["depth list"], D_confirmed)[5]
        model_fit_stats[8,i,j] = ksstats(summary_stat.stat["radius_ratio_list"].^2, D_ratio_confirmed)[5] #simulations save radius ratios while we computed transit duration ratios from the Kepler catalog

        println("# sigma_incl, sigma_hk: ", get_real(sim_param_here,"sigma_incl"), ", ",  get_real(sim_param_here,"sigma_hk"))
        toc()
    end
end
=#

#=
f = open("Model_stats_sigma_incl_sigma_hk_grid.out", "w")
println(f,"# Model statistics on a grid of (sigma_incl, sigma_hk) = (", sigma_incl_range, ", ",  sigma_hk_range, ")")
println(f, "# num_targets_sim_pass_one: ", get_int(sim_param,"num_targets_sim_pass_one"))
println(f, "# max_incl_sys: ", get_real(sim_param,"max_incl_sys"))
println(f, "# log_rate_clusters: ", get_real(sim_param,"log_rate_clusters"))
println(f, "# max_clusters_in_sys: ", get_int(sim_param,"max_clusters_in_sys"))
println(f, "# log_rate_planets_per_cluster: ", get_real(sim_param,"log_rate_planets_per_cluster"))
println(f, "# max_planets_in_clusters: ", get_int(sim_param,"max_planets_in_cluster"))
println(f, "# power_law_P: ", get_real(sim_param,"power_law_P"))
println(f, "# min_period: ", get_real(sim_param,"min_period"))
println(f, "# max_period: ", get_real(sim_param,"max_period"))
println(f, "# power_law_r1: ", get_real(sim_param,"power_law_r1"))
println(f, "# power_law_r2: ", get_real(sim_param,"power_law_r2"))
println(f, "# min_radius (R_earth): ", get_real(sim_param,"min_radius")/ExoplanetsSysSim.earth_radius)
println(f, "# max_radius (R_earth): ", get_real(sim_param,"max_radius")/ExoplanetsSysSim.earth_radius)
println(f, "# break_radius (R_earth): ", get_real(sim_param,"break_radius")/ExoplanetsSysSim.earth_radius)
#####println(f, "# sigma_incl: ", get_real(sim_param,"sigma_incl"))
println(f, "# sigma_incl_near_mmr: ", get_real(sim_param,"sigma_incl_near_mmr"))
#####println(f, "# sigma_hk: ", get_real(sim_param,"sigma_hk"))
println(f, "# num_mutual_hill_radii: ", get_real(sim_param,"num_mutual_hill_radii"))
println(f, "# mr_power_index: ", get_real(sim_param,"mr_power_index"))
println(f, "# mr_max_mass (M_earth): ", get_real(sim_param,"mr_max_mass")/ExoplanetsSysSim.earth_mass)
println(f, "# sigma_log_radius_in_cluster: ", get_real(sim_param,"sigma_log_radius_in_cluster"))
println(f, "# sigma_logperiod_per_pl_in_cluster: ", get_real(sim_param,"sigma_logperiod_per_pl_in_cluster"))
for i in 1:size(model_fit_stats)[1]
    println(f,"# ", model_fit_stats_fields[i], ":")
    println(f,model_fit_stats[i,:,:]')
end
close(f)
=#

#=
#To run our model on a 2D grid of 'mr_power_index' and 'sigma_log_radius_in_cluster':

mr_power_index_range = linspace(1.5, 3.5, 11)
sigma_log_radius_in_cluster_range = linspace(0.1, 1.0, 11)

model_fit_stats = zeros(8,length(mr_power_index_range),length(sigma_log_radius_in_cluster_range)) #3D array to be filled with all the statistics of the model fits: number of observed planets, and KS distances for multiplicities, periods, period ratios, transit durations, period-normalized transit duration ratios, transit depths, and transit depth ratios
model_fit_stats_fields = ["Number of observed planets", "KS multiplicities", "KS periods", "KS period ratios", "KS transit durations", "KS period-normalized transit duration ratios", "KS transit depths", "KS transit depth ratios"] #names of the fields associated with 'model_fit_stats' array

for (i,max_radius) in enumerate(mr_power_index_range)
    for (j,mr_power_index) in enumerate(sigma_log_radius_in_cluster_range)
        tic()
        sim_param_here = deepcopy(sim_param) #To set up all the default parameters
        add_param_fixed(sim_param_here,"mr_power_index",mr_power_index) #To replace the 'mr_power_index' parameter
        add_param_fixed(sim_param_here,"sigma_log_radius_in_cluster",sigma_log_radius_in_cluster) #To replace the 'sigma_log_radius_in_cluster' parameter

        cat_phys = generate_kepler_physical_catalog(sim_param_here)
        cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys,sim_param_here)
        cat_obs = observe_kepler_targets_single_obs(cat_phys_cut,sim_param_here)
        summary_stat = calc_summary_stats_model(cat_obs,sim_param_here)

        M_cat_obs = ones(Int64,0) #array to be filled with the number of transiting planets in each simulated system
        for k in 1:get_int(sim_param_here,"max_tranets_in_sys")
            append!(M_cat_obs, k*ones(Int64, summary_stat.stat["num n-tranet systems"][k]))
        end

        model_fit_stats[1,i,j] = summary_stat.stat["num_tranets"]
        model_fit_stats[2,i,j] = ksstats_ints(M_cat_obs, M_confirmed)[5]
        model_fit_stats[3,i,j] = ksstats(summary_stat.stat["P list"], P_confirmed)[5]
        model_fit_stats[4,i,j] = ksstats(summary_stat.stat["period_ratio_list"], R_confirmed)[5]
        model_fit_stats[5,i,j] = ksstats(summary_stat.stat["duration list"].*24, t_D_confirmed)[5] #transit durations in simulations are in days, while in the Kepler catalog are in hours
        model_fit_stats[6,i,j] = ksstats(summary_stat.stat["duration_ratio_list"], xi_confirmed)[5]
        model_fit_stats[7,i,j] = ksstats(summary_stat.stat["depth list"], D_confirmed)[5]
        model_fit_stats[8,i,j] = ksstats(summary_stat.stat["radius_ratio_list"].^2, D_ratio_confirmed)[5] #simulations save radius ratios while we computed transit duration ratios from the Kepler catalog

        println("# mr_power_index, sigma_log_radius_in_cluster: ", get_real(sim_param_here,"mr_power_index"), ", ",  get_real(sim_param_here,"sigma_log_radius_in_cluster"))
        toc()
    end
end
=#

#=
f = open("Model_stats_mr_sigmaR_grid.out", "w")
println(f,"# Model statistics on a grid of (mr_power_index,sigma_log_radius_in_cluster) = (", mr_power_index_range, ", ",  sigma_log_radius_in_cluster_range, ")")
println(f, "# num_targets_sim_pass_one: ", get_int(sim_param,"num_targets_sim_pass_one"))
println(f, "# max_incl_sys: ", get_real(sim_param,"max_incl_sys"))
println(f, "# log_rate_clusters: ", get_real(sim_param,"log_rate_clusters"))
println(f, "# max_clusters_in_sys: ", get_int(sim_param,"max_clusters_in_sys"))
println(f, "# log_rate_planets_per_cluster: ", get_real(sim_param,"log_rate_planets_per_cluster"))
println(f, "# max_planets_in_clusters: ", get_int(sim_param,"max_planets_in_cluster"))
println(f, "# power_law_P: ", get_real(sim_param,"power_law_P"))
println(f, "# min_period: ", get_real(sim_param,"min_period"))
println(f, "# max_period: ", get_real(sim_param,"max_period"))
println(f, "# power_law_r1: ", get_real(sim_param,"power_law_r1"))
println(f, "# power_law_r2: ", get_real(sim_param,"power_law_r2"))
println(f, "# min_radius (R_earth): ", get_real(sim_param,"min_radius")/ExoplanetsSysSim.earth_radius)
println(f, "# max_radius (R_earth): ", get_real(sim_param,"max_radius")/ExoplanetsSysSim.earth_radius)
println(f, "# break_radius (R_earth): ", get_real(sim_param,"break_radius")/ExoplanetsSysSim.earth_radius)
println(f, "# sigma_incl: ", get_real(sim_param,"sigma_incl"))
println(f, "# sigma_incl_near_mmr: ", get_real(sim_param,"sigma_incl_near_mmr"))
println(f, "# sigma_hk: ", get_real(sim_param,"sigma_hk"))
println(f, "# num_mutual_hill_radii: ", get_real(sim_param,"num_mutual_hill_radii"))
#####println(f, "# mr_power_index: ", get_real(sim_param,"mr_power_index"))
println(f, "# mr_max_mass (M_earth): ", get_real(sim_param,"mr_max_mass")/ExoplanetsSysSim.earth_mass)
#####println(f, "# sigma_log_radius_in_cluster: ", get_real(sim_param,"sigma_log_radius_in_cluster"))
println(f, "# sigma_logperiod_per_pl_in_cluster: ", get_real(sim_param,"sigma_logperiod_per_pl_in_cluster"))
for i in 1:size(model_fit_stats)[1]
    println(f,"# ", model_fit_stats_fields[i], ":")
    println(f,model_fit_stats[i,:,:]')
end
close(f)
=#

#=
#To run our model on a 2D grid of 'log_rate_clusters' and 'log_rate_planets_per_cluster':

log_rate_clusters_range = log.(linspace(1., 5., 11))
log_rate_planets_per_cluster_range = log.(linspace(1., 5., 11))

model_fit_stats = zeros(9,length(log_rate_clusters_range),length(log_rate_planets_per_cluster_range)) #3D array to be filled with all the statistics of the model fits: number of observed planets, and KS distances for multiplicities, periods, period ratios, transit durations, period-normalized transit duration ratios, transit depths, and transit depth ratios
model_fit_stats_fields = ["Number of observed planets", "Number of systems with observed singles", "KS multiplicities", "KS periods", "KS period ratios", "KS transit durations", "KS period-normalized transit duration ratios", "KS transit depths", "KS transit depth ratios"] #names of the fields associated with 'model_fit_stats' array

for (i,log_rate_clusters) in enumerate(log_rate_clusters_range)
    for (j,log_rate_planets_per_cluster) in enumerate(log_rate_planets_per_cluster_range)
        tic()
        sim_param_here = deepcopy(sim_param) #To set up all the default parameters
        add_param_fixed(sim_param_here,"log_rate_clusters",log_rate_clusters) #To replace the 'log_rate_clusters' parameter
        add_param_fixed(sim_param_here,"log_rate_planets_per_cluster",log_rate_planets_per_cluster) #To replace the 'log_rate_planets_per_cluster' parameter

        cat_phys = generate_kepler_physical_catalog(sim_param_here)
        cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys,sim_param_here)
        cat_obs = observe_kepler_targets_single_obs(cat_phys_cut,sim_param_here)
        summary_stat = calc_summary_stats_model(cat_obs,sim_param_here)

        M_cat_obs = ones(Int64,0) #array to be filled with the number of transiting planets in each simulated system
        for k in 1:get_int(sim_param_here,"max_tranets_in_sys")
            append!(M_cat_obs, k*ones(Int64, summary_stat.stat["num n-tranet systems"][k]))
        end

        model_fit_stats[1,i,j] = summary_stat.stat["num_tranets"]
        model_fit_stats[2,i,j] = summary_stat.stat["num n-tranet systems"][1]
        model_fit_stats[3,i,j] = ksstats_ints(M_cat_obs, M_confirmed)[5]
        model_fit_stats[4,i,j] = ksstats(summary_stat.stat["P list"], P_confirmed)[5]
        model_fit_stats[5,i,j] = ksstats(summary_stat.stat["period_ratio_list"], R_confirmed)[5]
        model_fit_stats[6,i,j] = ksstats(summary_stat.stat["duration list"].*24, t_D_confirmed)[5] #transit durations in simulations are in days, while in the Kepler catalog are in hours
        model_fit_stats[7,i,j] = ksstats(summary_stat.stat["duration_ratio_list"], xi_confirmed)[5]
        model_fit_stats[8,i,j] = ksstats(summary_stat.stat["depth list"], D_confirmed)[5]
        model_fit_stats[9,i,j] = ksstats(summary_stat.stat["radius_ratio_list"].^2, D_ratio_confirmed)[5] #simulations save radius ratios while we computed transit duration ratios from the Kepler catalog

        println("# exp(log_rate_clusters), exp(log_rate_planets_per_cluster): ", exp(get_real(sim_param_here,"log_rate_clusters")), ", ",  exp(get_real(sim_param_here,"log_rate_planets_per_cluster")))
        toc()
    end
end
=#

#=
f = open("Model_stats_Nc_Np_grid.out", "w")
println(f,"# Model statistics on a grid of (log_rate_clusters, log_rate_planets_per_cluster) = (", log_rate_clusters_range, ", ",  log_rate_planets_per_cluster_range, ")")
println(f, "# num_targets_sim_pass_one: ", get_int(sim_param,"num_targets_sim_pass_one"))
println(f, "# max_incl_sys: ", get_real(sim_param,"max_incl_sys"))
#####println(f, "# log_rate_clusters: ", get_real(sim_param,"log_rate_clusters"))
println(f, "# max_clusters_in_sys: ", get_int(sim_param,"max_clusters_in_sys"))
#####println(f, "# log_rate_planets_per_cluster: ", get_real(sim_param,"log_rate_planets_per_cluster"))
println(f, "# max_planets_in_clusters: ", get_int(sim_param,"max_planets_in_cluster"))
println(f, "# power_law_P: ", get_real(sim_param,"power_law_P"))
println(f, "# min_period: ", get_real(sim_param,"min_period"))
println(f, "# max_period: ", get_real(sim_param,"max_period"))
println(f, "# power_law_r1: ", get_real(sim_param,"power_law_r1"))
println(f, "# power_law_r2: ", get_real(sim_param,"power_law_r2"))
println(f, "# min_radius (R_earth): ", get_real(sim_param,"min_radius")/ExoplanetsSysSim.earth_radius)
println(f, "# max_radius (R_earth): ", get_real(sim_param,"max_radius")/ExoplanetsSysSim.earth_radius)
println(f, "# break_radius (R_earth): ", get_real(sim_param,"break_radius")/ExoplanetsSysSim.earth_radius)
println(f, "# sigma_incl: ", get_real(sim_param,"sigma_incl"))
println(f, "# sigma_incl_near_mmr: ", get_real(sim_param,"sigma_incl_near_mmr"))
println(f, "# sigma_hk: ", get_real(sim_param,"sigma_hk"))
println(f, "# num_mutual_hill_radii: ", get_real(sim_param,"num_mutual_hill_radii"))
println(f, "# mr_power_index: ", get_real(sim_param,"mr_power_index"))
println(f, "# mr_max_mass (M_earth): ", get_real(sim_param,"mr_max_mass")/ExoplanetsSysSim.earth_mass)
println(f, "# sigma_log_radius_in_cluster: ", get_real(sim_param,"sigma_log_radius_in_cluster"))
println(f, "# sigma_logperiod_per_pl_in_cluster: ", get_real(sim_param,"sigma_logperiod_per_pl_in_cluster"))
for i in 1:size(model_fit_stats)[1]
    println(f,"# ", model_fit_stats_fields[i], ":")
    println(f,model_fit_stats[i,:,:]')
end
close(f)
=#

##### The following optimization code is for exploring our updated model with a broken power law for the planet radii distribution instead of a single power law:

#=
#To run our model on a 2D grid of 'power_law_r1' and 'power_law_r2':

power_law_r1_range = linspace(-6., 3., 11)
power_law_r2_range = linspace(-6., 3., 11)

model_fit_stats = zeros(8,length(power_law_r1_range),length(power_law_r2_range)) #3D array to be filled with all the statistics of the model fits: number of observed planets, and KS distances for multiplicities, periods, period ratios, transit durations, period-normalized transit duration ratios, transit depths, and transit depth ratios
model_fit_stats_fields = ["Number of observed planets", "KS multiplicities", "KS periods", "KS period ratios", "KS transit durations", "KS period-normalized transit duration ratios", "KS transit depths", "KS transit depth ratios"] #names of the fields associated with 'model_fit_stats' array

for (i,power_law_r1) in enumerate(power_law_r1_range)
    for (j,power_law_r2) in enumerate(power_law_r2_range)
        tic()
        sim_param_here = deepcopy(sim_param) #To set up all the default parameters
        add_param_fixed(sim_param_here,"power_law_r1",power_law_r1) #To replace the 'power_law_r1' parameter
        add_param_fixed(sim_param_here,"power_law_r2",power_law_r2) #To replace the 'power_law_r2' parameter

        cat_phys = generate_kepler_physical_catalog(sim_param_here)
        cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys,sim_param_here)
        cat_obs = observe_kepler_targets_single_obs(cat_phys_cut,sim_param_here)
        summary_stat = calc_summary_stats_model(cat_obs,sim_param_here)

        M_cat_obs = ones(Int64,0) #array to be filled with the number of transiting planets in each simulated system
        for k in 1:get_int(sim_param_here,"max_tranets_in_sys")
            append!(M_cat_obs, k*ones(Int64, summary_stat.stat["num n-tranet systems"][k]))
        end

        model_fit_stats[1,i,j] = summary_stat.stat["num_tranets"]
        model_fit_stats[2,i,j] = ksstats_ints(M_cat_obs, M_confirmed)[5]
        model_fit_stats[3,i,j] = ksstats(summary_stat.stat["P list"], P_confirmed)[5]
        model_fit_stats[4,i,j] = ksstats(summary_stat.stat["period_ratio_list"], R_confirmed)[5]
        model_fit_stats[5,i,j] = ksstats(summary_stat.stat["duration list"].*24, t_D_confirmed)[5] #transit durations in simulations are in days, while in the Kepler catalog are in hours
        model_fit_stats[6,i,j] = ksstats(summary_stat.stat["duration_ratio_list"], xi_confirmed)[5]
        model_fit_stats[7,i,j] = ksstats(summary_stat.stat["depth list"], D_confirmed)[5]
        model_fit_stats[8,i,j] = ksstats(summary_stat.stat["radius_ratio_list"].^2, D_ratio_confirmed)[5] #simulations save radius ratios while we computed transit duration ratios from the Kepler catalog

        println("# power_law_r1, power_law_r2: ", get_real(sim_param_here,"power_law_r1"), ", ",  get_real(sim_param_here,"power_law_r2"))
        toc()
    end
end
=#

#=
f = open("Model_stats_r1_r2_rbreak2_grid.out", "w")
println(f,"# Model statistics on a grid of (power_law_r1, power_law_r2) = (", power_law_r1_range, ", ",  power_law_r2_range, ")")
println(f, "# num_targets_sim_pass_one: ", get_int(sim_param,"num_targets_sim_pass_one"))
println(f, "# max_incl_sys: ", get_real(sim_param,"max_incl_sys"))
println(f, "# log_rate_clusters: ", get_real(sim_param,"log_rate_clusters"))
println(f, "# max_clusters_in_sys: ", get_int(sim_param,"max_clusters_in_sys"))
println(f, "# log_rate_planets_per_cluster: ", get_real(sim_param,"log_rate_planets_per_cluster"))
println(f, "# max_planets_in_clusters: ", get_int(sim_param,"max_planets_in_cluster"))
println(f, "# power_law_P: ", get_real(sim_param,"power_law_P"))
println(f, "# min_period: ", get_real(sim_param,"min_period"))
println(f, "# max_period: ", get_real(sim_param,"max_period"))
#####println(f, "# power_law_r1: ", get_real(sim_param,"power_law_r1"))
#####println(f, "# power_law_r2: ", get_real(sim_param,"power_law_r2"))
println(f, "# min_radius (R_earth): ", get_real(sim_param,"min_radius")/ExoplanetsSysSim.earth_radius)
println(f, "# max_radius (R_earth): ", get_real(sim_param,"max_radius")/ExoplanetsSysSim.earth_radius)
println(f, "# break_radius (R_earth): ", get_real(sim_param,"break_radius")/ExoplanetsSysSim.earth_radius)
println(f, "# sigma_incl: ", get_real(sim_param,"sigma_incl"))
println(f, "# sigma_incl_near_mmr: ", get_real(sim_param,"sigma_incl_near_mmr"))
println(f, "# sigma_hk: ", get_real(sim_param,"sigma_hk"))
println(f, "# num_mutual_hill_radii: ", get_real(sim_param,"num_mutual_hill_radii"))
println(f, "# mr_power_index: ", get_real(sim_param,"mr_power_index"))
println(f, "# mr_max_mass (M_earth): ", get_real(sim_param,"mr_max_mass")/ExoplanetsSysSim.earth_mass)
println(f, "# sigma_log_radius_in_cluster: ", get_real(sim_param,"sigma_log_radius_in_cluster"))
println(f, "# sigma_logperiod_per_pl_in_cluster: ", get_real(sim_param,"sigma_logperiod_per_pl_in_cluster"))
for i in 1:size(model_fit_stats)[1]
    println(f,"# ", model_fit_stats_fields[i], ":")
    println(f,model_fit_stats[i,:,:]')
end
close(f)
=#

##### The following optimization code is for exploring our updated model with clustered (lognormal) planet radii, but still a single power law:

#=
#To run our model on a 2D grid of 'power_law_r' and 'sigma_log_radius_in_cluster':

power_law_r_range = linspace(-3., -1.2, 11)
sigma_log_radius_in_cluster_range = linspace(0.05, 0.5, 11)

model_fit_stats = zeros(7,length(power_law_r_range),length(sigma_log_radius_in_cluster_range)) #3D array to be filled with all the statistics of the model fits: number of observed planets, number of systems with observed singles, and KS distances for multiplicities, periods, period ratios, transit depths, and period-normalized transit duration ratios
model_fit_stats_fields = ["Number of observed planets", "Number of systems with observed singles", "KS multiplicities", "KS periods", "KS period ratios", "KS transit depths", "KS period-normalized transit duration ratios"] #names of the fields associated with 'model_fit_stats' array

sim_param_default = setup_sim_param_model() #To set up all the default parameters

for (i,power_law_r) in enumerate(power_law_r_range)
    for (j,sigma_log_radius_in_cluster) in enumerate(sigma_log_radius_in_cluster_range)
        tic()
        sim_param_here = deepcopy(sim_param_default) #To set up all the default parameters
        add_param_fixed(sim_param_here,"power_law_r",power_law_r) #To replace the 'power_law_r' parameter
        add_param_fixed(sim_param_here,"sigma_log_radius_in_cluster",sigma_log_radius_in_cluster) #To replace the 'sigma_log_radius_in_cluster' parameter

        cat_phys = generate_kepler_physical_catalog(sim_param_here)
        cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys,sim_param_here)
        cat_obs = observe_kepler_targets_single_obs(cat_phys_cut,sim_param_here)
        summary_stat = calc_summary_stats_model(cat_obs,sim_param_here)

        M_cat_obs = ones(Int64,0) #array to be filled with the number of transiting planets in each simulated system
        for k in 1:get_int(sim_param_here,"max_tranets_in_sys")
            append!(M_cat_obs, k*ones(Int64, summary_stat.stat["num n-tranet systems"][k]))
        end

        model_fit_stats[1,i,j] = summary_stat.stat["num_tranets"]
        model_fit_stats[2,i,j] = summary_stat.stat["num n-tranet systems"][1]
        model_fit_stats[3,i,j] = ksstats_ints(M_cat_obs, M_confirmed)[5]
        model_fit_stats[4,i,j] = ksstats(summary_stat.stat["P list"], P_confirmed)[5]
        model_fit_stats[5,i,j] = ksstats(summary_stat.stat["period_ratio_list"], R_confirmed)[5]
        model_fit_stats[6,i,j] = ksstats(summary_stat.stat["depth list"], D_confirmed)[5]
        model_fit_stats[7,i,j] = ksstats(summary_stat.stat["duration_ratio_list"], xi_confirmed)[5]

        println("# power_law_r, sigma_log_radius_in_cluster: ", get_real(sim_param_here,"power_law_r"), ", ",  get_real(sim_param_here,"sigma_log_radius_in_cluster"))
        toc()
    end
end
=#

#=
f = open("Model_stats_r_sigma_r_grid.out", "w")
println(f,"# Model statistics on a grid of (power_law_r, sigma_log_radius_in_cluster) = (", power_law_r_range, ", ",  sigma_log_radius_in_cluster_range, ")")
println(f,"# num_targets_sim_pass_one: ", get_int(sim_param_default,"num_targets_sim_pass_one"))
println(f,"# log_rate_clusters: ", get_real(sim_param_default,"log_rate_clusters"))
println(f,"# log_rate_planets_per_cluster: ", get_real(sim_param_default,"log_rate_planets_per_cluster"))
println(f,"# power_law_P: ", get_real(sim_param_default,"power_law_P"))
println(f,"# max_radius: ", get_real(sim_param_default,"max_radius"))
println(f,"# sigma_incl: ", get_real(sim_param_default,"sigma_incl"))
println(f,"# sigma_incl_near_mmr: ", get_real(sim_param_default,"sigma_incl_near_mmr"))
println(f,"# sigma_hk: ", get_real(sim_param_default,"sigma_hk"))
println(f,"# num_mutual_hill_radii: ", get_real(sim_param_default,"num_mutual_hill_radii"))
println(f,"# mr_power_index: ", get_real(sim_param_default,"mr_power_index"))
println(f,"# mr_max_mass: ", get_real(sim_param_default,"mr_max_mass"))
println(f,"# sigma_logperiod_per_pl_in_cluster: ", get_real(sim_param_default,"sigma_logperiod_per_pl_in_cluster"))
for i in 1:size(model_fit_stats)[1]
    println(f,"# ", model_fit_stats_fields[i], ":")
    println(f,model_fit_stats[i,:,:]')
end
close(f)
=#


text:
# VLSV reader in Julia

include("vlsvvariables.jl")

"Velocity mesh information."
struct VMeshInfo
   "number of velocity blocks"
   vblocks::SVector{3, Int64}
   vblock_size::SVector{3, Int64}
   vmin::SVector{3, Float64}
   vmax::SVector{3, Float64}
   dv::SVector{3, Float64}
end

"Variable MetaVLSV from the vlsv footer."
struct VarInfo
   "unit of the variable as a string"
   unit::String
   "unit of the variable as a LaTeX-formatted string"
   unitLaTeX::LaTeXString
   "description of the variable as a LaTeX-formatted string"
   variableLaTeX::LaTeXString
   "conversion factor to SI units as a string"
   unitConversion::String
end

"VLSV meta data."
struct MetaVLSV
   name::String
   dir::String
   fid::IOStream
   footer::EzXML.Node
   variable::Vector{String}
   "unsorted cell IDs"
   cellid::Vector{UInt64}
   "ordered sequence index of raw cell IDs"
   cellindex::Vector{Int64}
   time::Float64
   maxamr::Int64
   hasvdf::Bool
   ncells::SVector{3, Int64}
   block_size::SVector{3, Int64}
   coordmin::SVector{3, Float64}
   coordmax::SVector{3, Float64}
   dcoord::SVector{3, Float64}
   species::Vector{String}
   meshes::Dict{String, VMeshInfo}
end

function Base.show(io::IO, meta::MetaVLSV)
   println(io, "File: ", meta.name)
   println(io, "Time: ", round(meta.time, digits=2))
   println(io, "Dimension: $(ndims(meta))")
   println(io, "Maximum AMR level: $(meta.maxamr)")
   println(io, "Contains VDF: $(meta.hasvdf)")
   println(io, "variables: ", meta.variable)
end

function Base.show(io::IO, s::VarInfo)
   println(io, "Variable in LaTeX: ", s.variableLaTeX)
   println(io, "Unit: ", s.unit)
   println(io, "Unit in LaTeX: ", s.unitLaTeX)
   println(io, "Unit conversion: ", s.unitConversion)
end

"Return the xml footer of vlsv."
function getfooter(fid::IOStream)
   # First 8 bytes indicate big-endian or else
   endian_offset = 8
   seek(fid, endian_offset)
   # Obtain the offset of the XML file
   offset = read(fid, UInt64)
   seek(fid, offset)
   footer = read(fid, String) |> parsexml |> root
end

"Return size and type information for the object."
function getObjInfo(footer, name, tag, attr)
   local arraysize, datasize, datatype, vectorsize, variable_offset
   isFound = false

   for var in findall("//$tag", footer)
      if var[attr] == String(name)
         arraysize = parse(Int, var["arraysize"])
         datasize = parse(Int, var["datasize"])
         datatype = var["datatype"]
         vectorsize = parse(Int, var["vectorsize"])
         variable_offset = parse(Int, nodecontent(var))
         isFound = true
         break
      end
   end

   !isFound && throw(ArgumentError("unknown variable $name"))

   if datatype == "float"
      T = datasize == 4 ? Float32 : Float64
   elseif datatype == "int"
      T = datasize == 4 ? Int32 : Int64
   elseif datatype == "uint"
      T = datasize == 4 ? UInt32 : UInt64
   end

   T, variable_offset, arraysize, datasize, vectorsize
end

"Return vectors of `name` from the vlsv file with `footer` opened by `fid`."
function readvector(fid::IOStream, footer, name, tag)
   T, offset, arraysize, datasize, vectorsize = getObjInfo(footer, name, tag, "name")

   if Sys.total_memory() > 8*arraysize*vectorsize*datasize
      w = vectorsize == 1 ?
         Vector{T}(undef, arraysize) :
         Array{T,2}(undef, vectorsize, arraysize)
      seek(fid, offset)
      read!(fid, w)
   else
      @warn "Large array detected. Using memory-mapped I/O!" maxlog=1
      a = mmap(fid, Vector{UInt8}, sizeof(T)*vectorsize*arraysize, offset)
      w = vectorsize == 1 ?
         reinterpret(T, a) :
         reshape(reinterpret(T, a), vectorsize, arraysize)
   end

   w
end

"""
    load(file) -> MetaVLSV

Return MetaVLSV from a vlsv `file`.
"""
function load(file::AbstractString)
   isfile(file) || throw(ArgumentError("Cannot open \'$file\': not a file"))
   fid = open(file, "r")

   footer = getfooter(fid)

   local cellid
   let
      T, offset, arraysize, _, vectorsize = getObjInfo(footer, "CellID", "VARIABLE", "name")
      a = mmap(fid, Vector{UInt8}, sizeof(T)*vectorsize*arraysize, offset)
      cellid = reinterpret(T, a)
   end

   cellindex = sortperm(cellid)

   bbox = readmesh(fid, footer, "SpatialGrid", "MESH_BBOX")::Vector{UInt}

   nodeCoordsX = readmesh(fid, footer, "SpatialGrid", "MESH_NODE_CRDS_X")::Vector{Float64}
   nodeCoordsY = readmesh(fid, footer, "SpatialGrid", "MESH_NODE_CRDS_Y")::Vector{Float64}
   nodeCoordsZ = readmesh(fid, footer, "SpatialGrid", "MESH_NODE_CRDS_Z")::Vector{Float64}

   @inbounds ncells = SVector(bbox[1], bbox[2], bbox[3])
   @inbounds block_size = SVector(bbox[4], bbox[5], bbox[6])
   @inbounds coordmin = SVector(nodeCoordsX[begin], nodeCoordsY[begin], nodeCoordsZ[begin])
   @inbounds coordmax = SVector(nodeCoordsX[end], nodeCoordsY[end], nodeCoordsZ[end])

   dcoord = SVector{3}(@. (coordmax - coordmin) / ncells)

   meshes = Dict{String, VMeshInfo}()

   # Find all species by the BLOCKIDS tag
   species = String[]

   for varinfo in findall("//BLOCKIDS", footer)

      if haskey(varinfo, "name")
         # VLSV 5.0 file with bounding box
         popname = varinfo["name"]

         bbox = readmesh(fid, footer, popname, "MESH_BBOX")::Vector{UInt}

         nodeCoordsX = readmesh(fid, footer, popname, "MESH_NODE_CRDS_X")::Vector{Float64}
         nodeCoordsY = readmesh(fid, footer, popname, "MESH_NODE_CRDS_Y")::Vector{Float64}
         nodeCoordsZ = readmesh(fid, footer, popname, "MESH_NODE_CRDS_Z")::Vector{Float64}
         vblocks = SVector(bbox[1], bbox[2], bbox[3])
         vblock_size = SVector(bbox[4], bbox[5], bbox[6])
         vmin = SVector(nodeCoordsX[begin], nodeCoordsY[begin], nodeCoordsZ[begin])
         vmax = SVector(nodeCoordsX[end], nodeCoordsY[end], nodeCoordsZ[end])
         dv = SVector{3}(@. (vmax - vmin) / vblocks / vblock_size)
      else
         popname = "avgs"

         if "vxblocks_ini" in getindex.(findall("//PARAMETER", footer), "name")
            # In VLSV before 5.0 the mesh is defined with parameters.
            vblocks = SVector(
               readparameter(fid, footer, "vxblocks_ini"),
               readparameter(fid, footer, "vyblocks_ini"),
               readparameter(fid, footer, "vzblocks_ini"))
            vblock_size = SVector(4, 4, 4)
            vmin = SVector(
               readparameter(fid, footer, "vxmin"),
               readparameter(fid, footer, "vymin"),
               readparameter(fid, footer, "vzmin") )
            vmax = SVector(
               readparameter(fid, footer, "vxmax"),
               readparameter(fid, footer, "vymax"),
               readparameter(fid, footer, "vzmax") )
            dv = SVector{3}(@. (vmax - vmin) / vblocks / vblock_size)
         else
            # No velocity space info, e.g., file not written by Vlasiator
            vblocks = SVector(0, 0, 0)
            vblock_size = SVector(4, 4, 4)
            vmin = SVector(0, 0, 0)
            vmax = SVector(0, 0, 0)
            dv = SVector(1, 1, 1)
         end
      end

      # Update list of active species
      if popname ∉ species
         push!(species, popname)
      end

      # Create a new object for this population
      popVMesh = VMeshInfo(vblocks, vblock_size, vmin, vmax, dv)

      meshes[popname] = popVMesh
   end

   if hasname(footer, "PARAMETER", "time") # Vlasiator 5.0+
      timesim = readparameter(fid, footer, "time")::Float64
   elseif hasname(footer, "PARAMETER", "t")
      timesim = readparameter(fid, footer, "t")::Float64
   else
      timesim = Inf
   end

   # Obtain maximum refinement level
   ncell = prod(ncells)
   maxamr, cid = 0, ncell
   while @inbounds cid < cellid[cellindex[end]]
      maxamr += 1
      cid += ncell*8^maxamr
   end

   varinfo = findall("//VARIABLE", footer)
   nVar = length(varinfo)
   vars = Vector{String}(undef, nVar)
   for i in eachindex(vars, varinfo)
      @inbounds vars[i] = varinfo[i]["name"]
   end

   hasvdf = let
      vcells = readmesh(fid, footer, "SpatialGrid", "CELLSWITHBLOCKS")::Vector{UInt}
      !isempty(vcells)
   end

   # File IOstream is not closed for sake of data processing later.

   meta = MetaVLSV(basename(file), dirname(file), fid, footer, vars, cellid,
      cellindex, timesim, maxamr, hasvdf, ncells, block_size, coordmin, coordmax,
      dcoord, species, meshes)
end

"""
    readvariablemeta(meta, var) -> VarInfo

Return VarInfo about `var` in the vlsv file linked to `meta`.
"""
function readvariablemeta(meta::MetaVLSV, var)

   varSym = isa(var, AbstractString) ? Symbol(var) : var

   unit, unitLaTeX, variableLaTeX, unitConversion = "", "", "", ""

   if varSym in keys(units_predefined)
      unit, variableLaTeX, unitLaTeX = units_predefined[varSym]
   elseif hasvariable(meta, var) # For Vlasiator 5 vlsv files, MetaVLSV is included
      for varinfo in findall("//VARIABLE", meta.footer)
         if varinfo["name"] == var
            unit = varinfo["unit"]
            unitLaTeX = varinfo["unitLaTeX"]
            variableLaTeX = varinfo["variableLaTeX"]
            unitConversion = varinfo["unitConversion"]
         end
         # If var isn't predefined or unit isn't found, it will return nothing!
      end
   end

   VarInfo(unit, unitLaTeX, variableLaTeX, unitConversion)
end

"Return mesh related variable."
function readmesh(fid::IOStream, footer, typeMesh, varMesh)
   T, offset, arraysize, _, _ = getObjInfo(footer, typeMesh, varMesh, "mesh")

   w = Vector{T}(undef, arraysize)
   seek(fid, offset)
   read!(fid, w)

   w
end

"""
    readvariable(meta::MetaVLSV, var, sorted::Bool=true) -> Array

Return variable value of `var` from the vlsv file. By default `sorted=true`, which means
that for DCCRG grid the variables are sorted by cell ID.
"""
function readvariable(meta::MetaVLSV, var, sorted::Bool=true)
   @unpack fid, footer, cellindex = meta
   if (local symvar = Symbol(var)) in keys(variables_predefined)
      data = variables_predefined[symvar](meta)
      return data
   end

   raw = readvector(fid, footer, var, "VARIABLE")

   if startswith(var, "fg_") # fsgrid
      bbox = SVector{6}(readmesh(fid, footer, "fsgrid", "MESH_BBOX")::Vector{Int})
      # Determine fsgrid domain decomposition
      nIORanks = readparameter(meta, "numWritingRanks")::Int32

      if ndims(raw) > 1
         @inbounds dataOrdered = zeros(Float32, size(raw,1), bbox[1], bbox[2], bbox[3])
      else
         @inbounds dataOrdered = zeros(Float32, bbox[1], bbox[2], bbox[3])
      end

      @inbounds fgDecomposition = @views getDomainDecomposition(bbox[1:3], nIORanks)

      currentOffset = ones(UInt32, nIORanks+1)
      lsize = ones(Int, 3, nIORanks)
      lstart = similar(lsize)
      @inbounds @views @floop for i = 1:nIORanks
         xyz = SA[
            (i - 1) ÷ fgDecomposition[3] ÷ fgDecomposition[2],
            (i - 1) ÷ fgDecomposition[3] % fgDecomposition[2],
            (i - 1) % fgDecomposition[3] ]

         lsize[:,i] = calcLocalSize.(bbox[1:3], fgDecomposition, xyz)
         lstart[:,i] = calcLocalStart.(bbox[1:3], fgDecomposition, xyz)

         totalSize = prod(lsize[:,i])
         currentOffset[i+1] = currentOffset[i] + totalSize

         # Reorder data
         if ndims(raw) > 1
            lend = lstart[:,i] + lsize[:,i] .- 1

            ldata = raw[:,currentOffset[i]:currentOffset[i+1]-1]
            ldata = reshape(ldata, size(raw,1), lsize[1,i], lsize[2,i], lsize[3,i])

            dataOrdered[:,lstart[1,i]:lend[1],lstart[2,i]:lend[2],lstart[3,i]:lend[3]] =
               ldata
         else
            ldata = raw[currentOffset[i]:currentOffset[i+1]-1]
            ldata = reshape(ldata, lsize[1,i], lsize[2,i], lsize[3,i])

            dataOrdered[lstart[1,i]:lend[1],lstart[2,i]:lend[2],lstart[3,i]:lend[3]] = ldata
         end
      end

      data = dropdims(dataOrdered, dims=(findall(size(dataOrdered) .== 1)...,))
   elseif sorted # dccrg grid
      @inbounds d = ndims(raw) == 1 ? raw[cellindex] : raw[:,cellindex]
      data = eltype(raw) == Float64 ? Float32.(d) : d
   else
      data = raw
   end
   return data
end

"""
    readvariable(meta::MetaVLSV, var, ids) -> Array

Read a variable `var` in a collection of cells `ids`.
"""
function readvariable(meta::MetaVLSV, var, ids)
   @assert !startswith(var, "fg_") "Currently does not support reading fsgrid!"
   @unpack fid, footer, cellid, cellindex = meta

   if (local symvar = Symbol(var)) in keys(variables_predefined)
      data = variables_predefined[symvar](meta, ids)
      return data
   end

   T, offset, arraysize, _, vectorsize = getObjInfo(footer, var, "VARIABLE", "name")

   v = Array{T}(undef, vectorsize, length(ids))

   a = mmap(fid, Vector{UInt8}, sizeof(T)*vectorsize*arraysize, offset)
   w = reshape(reinterpret(T, a), vectorsize, arraysize)

   id_ = length(cellindex) < 1000 ?
      [findfirst(==(id), cellid) for id in ids] :
      indexin(ids, cellid)

   for i in eachindex(id_), iv = 1:vectorsize
      @inbounds v[iv,i] = w[iv,id_[i]]
   end
   if T == Float64
      v = Float32.(v)
   end
   return v
end

@inline Base.getindex(meta::MetaVLSV, key::AbstractString) = readvariable(meta, key)

"Return 2d scalar/vector data. Nonpublic because it won't work with DCCRG AMR."
function getdata2d(meta::MetaVLSV, var)
   @assert ndims(meta) == 2 "2D outputs required."
   sizes = filter(!=(1), meta.ncells)
   data = readvariable(meta, var)
   data = ndims(data) == 1 ?
      reshape(data, sizes[1], sizes[2]) :
      reshape(data, 3, sizes[1], sizes[2]) # assumes 3-vector, may not work in general
end

"File size in bytes."
@inline Base.size(meta::MetaVLSV) = filesize(joinpath(meta.dir, meta.name))

# Optimize decomposition of this grid over the given number of processors.
# Reference: fsgrid.hpp
function getDomainDecomposition(globalsize, nprocs)
   domainDecomp = SVector(1, 1, 1)
   minValue = typemax(Int)

   @inbounds for i = 1:min(nprocs, globalsize[1])
      iBox = max(globalsize[1]/i, 1)

      for j = 1:min(nprocs, globalsize[2])
         i * j > nprocs && break

         jBox = max(globalsize[2]/j, 1)

         for k = 1:min(nprocs, globalsize[2])
            i * j * k > nprocs && continue

            kBox = max(globalsize[3]/k, 1)

            nyz = ifelse(i > 1, jBox * kBox, 0)
            nzx = ifelse(j > 1, kBox * iBox, 0)
            nxy = ifelse(k > 1, iBox * jBox, 0)

            v = 10*iBox*jBox*kBox + nyz + nzx + nxy

            if i * j * k == nprocs && v < minValue
               minValue = v
               domainDecomp = SVector(i, j, k)
            end
         end
      end
   end

   domainDecomp
end

function calcLocalStart(globalCells, nprocs, lcells)
   ncells = globalCells ÷ nprocs
   remainder = globalCells % nprocs
   lstart = ifelse(lcells < remainder, lcells*(ncells+1) + 1, lcells*ncells + remainder + 1)
end

function calcLocalSize(globalCells, nprocs, lcells)
   ncells = globalCells ÷ nprocs
   remainder = globalCells % nprocs
   lsize = ifelse(lcells < remainder, ncells + 1, ncells)
end

"""
    hasvariable(meta, var) -> Bool

Check if the VLSV file contains a variable.
"""
hasvariable(meta::MetaVLSV, var) = hasname(meta.footer, "VARIABLE", var)

"""
    readparameter(meta, param)

Return the parameter value from vlsv file.
"""
readparameter(meta::MetaVLSV, param) = readparameter(meta.fid, meta.footer, param)

function readparameter(fid::IOStream, footer, param)
   T, offset, _, _, _ = getObjInfo(footer, param, "PARAMETER", "name")
   seek(fid, offset)
   p = read(fid, T)
end

"""
    hasparameter(meta, param) -> Bool

Check if the vlsv file contains a certain parameter.
"""
hasparameter(meta::MetaVLSV, param) = hasname(meta.footer, "PARAMETER", param)

"Check if the XML `element` contains a `tag` with `name`."
function hasname(element, tag, name)
   isFound = false

   for var in findall("//$tag", element)
      var["name"] == name && (isFound = true)
      isFound && break
   end

   isFound
end

"""
    ndims(meta) -> Int

Return the dimension of VLSV data.
"""
Base.ndims(meta::MetaVLSV) = count(>(1), meta.ncells)

"""
    readvcells(meta, cid; species="proton")

Read velocity cells from a spatial cell of ID `cid`, and return a map of velocity cell
ids and corresponding value.
"""
function readvcells(meta::MetaVLSV, cid; species="proton")
   @unpack fid, footer = meta
   @unpack vblock_size = meta.meshes[species]
   bsize = prod(vblock_size)

   local offset, nblocks
   let cellsWithVDF = readvector(fid, footer, species, "CELLSWITHBLOCKS"),
       nblock_C = readvector(fid, footer, species, "BLOCKSPERCELL")
      # Check if cells have vspace stored
      if cid ∈ cellsWithVDF
         cellWithVDFIndex = findfirst(==(cid), cellsWithVDF)
      else
         throw(ArgumentError("Cell ID $cid does not store velocity distribution!"))
      end
      # Navigate to the correct position
      offset = sum(@view nblock_C[1:cellWithVDFIndex-1])
      @inbounds nblocks = Int(nblock_C[cellWithVDFIndex])
   end

   local datasize, datatype, vectorsize, variable_offset
   # Read in avgs
   let varinfo = findfirst("//BLOCKVARIABLE", footer)
      datasize = parse(Int, varinfo["datasize"])
      datatype = varinfo["datatype"]
      datatype == "float" || @error "VDFs must be floating numbers!"
      vectorsize = parse(Int, varinfo["vectorsize"])
      variable_offset = parse(Int, nodecontent(varinfo))
   end

   data = let
      Tavg = datasize == 4 ? Float32 : Float64
      a = mmap(fid, Vector{UInt8}, sizeof(Tavg)*vectorsize*nblocks,
         offset*vectorsize*datasize + variable_offset)
      reshape(reinterpret(Tavg, a), vectorsize, nblocks)
   end

   # Read in block IDs
   let varinfo = findfirst("//BLOCKIDS", footer)
      datasize = parse(Int, varinfo["datasize"])
      datatype = varinfo["datatype"]
      datatype == "uint" || @error "Block ID must be unsigned integer!"
      variable_offset = parse(Int, nodecontent(varinfo))
   end

   blockIDs = let T = datasize == 4 ? UInt32 : UInt64
      Vector{T}(undef, nblocks)
   end
   seek(fid, offset * datasize + variable_offset)
   read!(fid, blockIDs)

   # Velocity cell IDs and distributions (ordered by blocks)
   vcellids = zeros(UInt32, bsize*nblocks)
   vcellf = zeros(Float32, bsize*nblocks)

   vcellid_local = @inbounds [i + vblock_size[1]*j + vblock_size[1]*vblock_size[2]*k
      for i in 0:vblock_size[1]-1, j in 0:vblock_size[2]-1, k in 0:vblock_size[3]-1]

   @inbounds @floop for i in eachindex(blockIDs), j = 1:bsize
      vcellids[(i-1)*bsize+j] = vcellid_local[j] + bsize*blockIDs[i]
      vcellf[(i-1)*bsize+j] = data[j,i]
   end
   vcellids, vcellf
end

ext: jl


text:
function convert(obj::PropagatorConverter, arg0::List, arg1::jboolean, arg2::List)
    return jcall(obj, "convert", Propagator, (List, jboolean, List), arg0, arg1, arg2)
end

function convert(obj::PropagatorConverter, arg0::List, arg1::jboolean, arg2::Vector{JString})
    return jcall(obj, "convert", Propagator, (List, jboolean, Vector{JString}), arg0, arg1, arg2)
end

function convert(obj::PropagatorConverter, arg0::Propagator, arg1::jdouble, arg2::jint, arg3::List)
    return jcall(obj, "convert", Propagator, (Propagator, jdouble, jint, List), arg0, arg1, arg2, arg3)
end

function convert(obj::PropagatorConverter, arg0::Propagator, arg1::jdouble, arg2::jint, arg3::Vector{JString})
    return jcall(obj, "convert", Propagator, (Propagator, jdouble, jint, Vector{JString}), arg0, arg1, arg2, arg3)
end


text:
module MDToolbox

using FileIO
using NetCDF
#using FFTW

# satndard library
using Printf
using Statistics
using LinearAlgebra
using Base.Threads

# package code goes here
export KB_kcalpermol, KB_kjpermol
export AbstractTrajectory, TrjArray, select_atom
export readdcd, readnetcdf, writenetcdf, readpsf, writepsf, readpdb, writepdb
export centerofmass, decenter, superimpose, superimpose_serial, getrmsd, meanstructure, getrmsf
export getdistance, getangle, getdihedral
export ksdensity, ksdensity_serial, getpmf
export wham

# constants
const KB_kcalpermol = 0.0019872041 #Boltzmann constant taken from wikipedia
const KB_kjpermol = 0.0083144621 #Boltzmann constant taken from wikipedia

# codes
include("trjarray.jl")
include("fileio.jl")
include("structure.jl")
include("ksdensity.jl")
include("wham.jl")

end # module


text:
include("helpers.jl")
include("Blocks/Sources.jl")

import ControlSystems

# Alias to make extensions less verbose. Maybe I should import the functions?
CS = ControlSystems

# ---------------------------- Generalized state space--------------------------------------
# Generalized state space form:
#       ẋ = f(x, u, t)
#       y = h(x, u, t)
#
#   where
#       x: internal states variables
#       u: input state variables
#       y: output state variables
#       t: simulation time
#       f: state function
#       h: output function
# ------------------------------------------------------------------------------------------
abstract type SSFunction end

# State function (f from generalized state space form)
struct StateFunction{Initialized} <: SSFunction
    func::Function
    init_cond

    StateFunction(func, init_cond) = new{true}(func, init_cond)
    StateFunction(func) = new{false}(func, nothing)
end

# Output function (h from generalized state space form)
struct OutputFunction{Feedthrough} <: SSFunction
    func::Function

    function OutputFunction(func)
        narg_vect = nargin(func)
        nmethods = length(narg_vect)
        nargs = narg_vect[1]

        @assert(nmethods==1,
            "Output function must have only 1 method, instead it has $nmethods"
        )

        if nargs==2
            feedthrough = false
        elseif nargs==3
            feedthrough = true
        else
            error("Output function must have 2 or 3 arguments, instead it has $nargs")
        end

        return new{feedthrough}(func)
    end
end

# # Feedthrough output functions are of form:
# #   h(x, u, t)
# struct FeedthroughOutputFunction <: AbstractOutputFunction
#     func::Function
#     nu::Int
#     ny::Int
# end
#
# # Output functions (non-feedthrough) are of form:
# #   h(x, t)
# struct OutputFunction <: AbstractOutputFunction
#     func::Function
#     ny::Int
# end

(f::SSFunction)(args...; kwargs...) = f.func(args...; kwargs...)

# function CS.state_space_validation(f::Function, h::Function, Ts)
#     f_narg_vect, h_narg_vect = (f, h) .|> nargin
#
#     @assert(length(f_narg_vect)==length(h_narg_vect)==1,
#         "State and output functions must have only one method"
#     )
#     @assert(Ts ≥ 0 || Ts == -1,
#         "Ts must be either a positive number, 0 (continuous system), or -1 (unspecified)"
#     )
#
#     f_nargs, h_nargs = f_narg_vect[1], h_narg_vect[1]
#
#
#     if h_nargs==2
#         h_type = OutputFunction
#     elseif h_nargs==3
#         h_type = FeedthroughOutputFunction
#     else
#         error("Output function must have signature: h(x, u, t) or h(x, t)")
#     end
#
#     return f.nx, f_nargs, h_nargs, h_type
# end

struct NLStateSpace{Initialized, Feedthrough}
    f::StateFunction{Initialized}
    h::OutputFunction{Feedthrough}
    x
    y
    Ts
    in
    out
end
function NLStateSpace(f, h, Ts=-1, in=[], out=[])
    f = StateFunction(f)
    h = OutputFunction(h)
    return NLStateSpace(f, h, Ts, in, out)
end

# abstract type AbstractNLStateSpace end
#
# struct FeedthroughNLStateSpace <: AbstractNLStateSpace
#     f::StateFunction
#     h::FeedthroughOutputFunction
#     Ts #::Float64
#     in
#     out
# end
#
# struct NLStateSpace <: AbstractNLStateSpace
#     f::StateFunction
#     h::OutputFunction
#     Ts #::Float64
#     in
#     out
#     NLStateSpace(f, h, Ts=-1) = new(f, h, Ts, CS.state_space_validation(f, h, Ts)...)
# end
#
#
# function (sys::AbstractNLStateSpace)(u)
#
# end
#
# # Extend ss method to create NLStateSpace type when functions are passed in
# function CS.ss(f::Function, h::Function, args..., kwargs...)
#     return
# end
# CS.ss(f::Function, h::OutputFunction, args...; kwargs...) =
#     NLStateSpace(f, h, args...; kwargs...)
# CS.ss(f::Function, h::FeedthroughOutputFunction, args...; kwargs...) =
#     FeedthroughNLStateSpace(f, h, args...; kwargs...)
#
# GeneralizedStateSpace = Union{CS.StateSpace, NLStateSpace}

initial_conditions(f::StateFunction) = f.init_cond
initial_conditions(sys::NLStateSpace) = sys.f.init_cond

CS.ninputs(sys::NLStateSpace) = sys.nu
CS.noutputs(sys::NLStateSpace) = sys.ny

# Indexing Functions
Base.ndims(::NLStateSpace) = 2
Base.size(sys::NLStateSpace) = (noutputs(sys), ninputs(sys))
Base.size(sys::NLStateSpace, d::Integer) = d <= 2 ? size(sys)[d] : 1
Base.eltype(::Type{S}) where {S<:NLStateSpace} = S

max_stars_repo_name: jonniedie/ComposableDiffEqs

max_issues_repo_name: jonniedie/ComposableDiffEqs

max_forks_repo_name: jonniedie/ComposableDiffEqs


text:
#=

# 202: 2D Laplace equation on unstructured mesh
([source code](SOURCE_URL))

=#
module Example202_Laplace2D_Unstructured

using VoronoiFVM
using LinearAlgebra

# Flux function which describes the flux
# between neigboring control volumes $\omega_k$ and $\omega_l$
function flux!(f,u,edge)
    f[1]=u[1,1]-u[1,2]
end

function source!(f,edge)
    f[1]=1
end

function main(;Plotter=nothing, plot_grid=false,nref=0)

    nspecies=1
    ispec=1
    nrad=10*2^nref
    grid=VoronoiFVM.Grid(points=reduce(hcat,[ [cos(2*pi*i/nrad), sin(2*pi*i/nrad)] for i=1:nrad]),
                     bfaces=reduce(hcat,vcat([ [i,i+1] for i=1:nrad-1],[[nrad,1]])),
                     bfaceregions=ones(nrad),
                     regionpoints=[0.0 0.0;],
                     regionnumbers=[1],
                     regionvolumes=[0.1*2.0^(-2*nref)])

    if plot_grid
        p=plot(Plotter=Plotter,grid)
        return
    end

    # Create a physics structure
    physics=VoronoiFVM.Physics(num_species=nspecies,flux=flux!,source=source!)

    # Create a finite volume system with dense storage of unknowns
    sys=VoronoiFVM.System(grid,physics,unknown_storage=:dense)

    # Enable species 1 in region 1
    enable_species!(sys,ispec,[1])

    # Set boundary conditions
    # Dirichlet boundary conditions are marked by setting a corresponding value of the
    # boundary factor
    for i=1:num_bfaceregions(grid)
        boundary_dirichlet!(sys,ispec,i,0.0)
    end

    # Create & initialize array for solution and initial value
    inival=unknowns(sys)
    solution=unknowns(sys)
    inival.=0
    solution.=0

    # Solve stationary problem
    solve!(solution,inival,sys)

    if Plotter!=nothing
        plot(grid, solution[1,:], Plotter=Plotter)
    end

    # Return test value
    return norm(solution,Inf)
end

#
# Called by unit test
#
function test()
    main() ≈ 0.23422190956113936
end

end


text:
using LinearAlgebra, CUDA, Flux, FastSolversForWeightedTV, Test
CUDA.allowscalar(false)

# Geometry
T = Float32
n2d = (256, 256)
h2d = (abs(randn(T)), abs(randn(T)))
n3d = (256, 256, 256)
h3d = (abs(randn(T)), abs(randn(T)), abs(randn(T)))
nc = 3
nb = 16

# Operators
flag_gpu = true
# flag_gpu = false
∇2d    = gradient_operator(n2d, h2d; T=T);                        flag_gpu && (∇2d    = ∇2d    |> gpu)
∇2d_c  = gradient_operator(n2d, h2d; T=Complex{T});               flag_gpu && (∇2d_c  = ∇2d_c  |> gpu)
∇2db   = gradient_operator_batch(n2d, nc, nb, h2d; T=T);          flag_gpu && (∇2db   = ∇2db   |> gpu)
∇2db_c = gradient_operator_batch(n2d, nc, nb, h2d; T=Complex{T}); flag_gpu && (∇2db_c = ∇2db_c |> gpu)
∇3d    = gradient_operator(n3d, h3d; T=T);                        flag_gpu && (∇3d    = ∇3d    |> gpu)
∇3d_c  = gradient_operator(n3d, h3d; T=Complex{T});               flag_gpu && (∇3d_c  = ∇3d_c  |> gpu)

# Adjoint test (∇2d)
u = randn(T, n2d);                   flag_gpu && (u = u |> gpu)
v = randn(T, n2d[1]-1, n2d[2]-1, 2); flag_gpu && (v = v |> gpu)
a = dot(∇2d*u, v)
b = dot(u, adjoint(∇2d)*v)
@test a ≈ b rtol = T(1e-5)

# Adjoint test (∇2d_c)
u = randn(Complex{T}, n2d);                   flag_gpu && (u = u |> gpu)
v = randn(Complex{T}, n2d[1]-1, n2d[2]-1, 2); flag_gpu && (v = v |> gpu)
a = dot(∇2d_c*u, v)
b = dot(u, adjoint(∇2d_c)*v)
@test a ≈ b rtol = T(1e-5)

# Adjoint test (∇2db)
u = randn(T, n2d..., nc, nb);               flag_gpu && (u = u |> gpu)
v = randn(T, n2d[1]-1, n2d[2]-1, 2*nc, nb); flag_gpu && (v = v |> gpu)
a = dot(∇2db*u, v)
b = dot(u, adjoint(∇2db)*v)
@test a ≈ b rtol = T(1e-5)

# Adjoint test (∇2db_c)
u = randn(Complex{T}, n2d..., nc, nb);               flag_gpu && (u = u |> gpu)
v = randn(Complex{T}, n2d[1]-1, n2d[2]-1, 2*nc, nb); flag_gpu && (v = v |> gpu)
a = dot(∇2db_c*u, v)
b = dot(u, adjoint(∇2db_c)*v)
@test a ≈ b rtol = T(1e-5)

# Adjoint test (∇3d)
u = randn(T, n3d);                             flag_gpu && (u = u |> gpu)
v = randn(T, n3d[1]-1, n3d[2]-1, n3d[3]-1, 3); flag_gpu && (v = v |> gpu)
a = dot(∇3d*u, v)
b = dot(u, adjoint(∇3d)*v)
@test a ≈ b rtol = T(1e-5)

# Adjoint test (∇3d_c)
u = randn(Complex{T}, n3d);                             flag_gpu && (u = u |> gpu)
v = randn(Complex{T}, n3d[1]-1, n3d[2]-1, n3d[3]-1, 3); flag_gpu && (v = v |> gpu)
a = dot(∇3d_c*u, v)
b = dot(u, adjoint(∇3d_c)*v)
@test a ≈ b rtol = T(1e-5)

max_stars_repo_name: grizzuti/FastSolversForWeightedTV

max_issues_repo_name: grizzuti/FastSolversForWeightedTV

max_forks_repo_name: grizzuti/FastSolversForWeightedTV


text:
include("../mpiops.jl")

mpi_init()

function kappa(x, y)
    return 1+x^2+y^2
end

function rhs(x, y)
    -2*x*(1 - x)*(x^2 + y^2 + 1) + 2*x*(-x*y*(1 - y) + y*(1 - x)*(1 - y)) -
             2*y*(1 - y)*(x^2 + y^2 + 1) + 2*y*(-x*y*(1 - x) + x*(1 - x)*(1 - y))
end

function ufunc(x, y)
    x * (1-x) * y * (1-y)
end

mc = MPIConfig(50)
global_to_local, local_to_global = dofmap(mc)
X, Y = get_xy(mc)
f_local = rhs.(X, Y)
using Random; Random.seed!(233)
θ = Variable(fc_init([2,20,20,20,1]))
θ = mpi_bcast(θ)
κ_local = (fc([X'[:] Y'[:]], [20,20,20,1], θ) + 5.0)|>squeeze
κ_local = reshape(κ_local, (mc.n, mc.n))
u_local = poisson_solver(κ_local, f_local, mc)
u = mpi_gather(u_local)[global_to_local]

Uexact = reshape(constant(ufunc.(X, Y)), (-1,))
Uexact = mpi_gather(Uexact)[global_to_local]

sess = Session(); init(sess)

U = run(sess, u)
Uexact = run(sess, Uexact)

if mpi_rank()==0
    close("all")
    figure(figsize = (15,4))
    subplot(131)
    pcolormesh(reshape(U, mc.n * mc.N, mc.n * mc.N))
    colorbar()
    subplot(132)
    pcolormesh(reshape(Uexact, mc.n * mc.N, mc.n * mc.N))
    colorbar()
    subplot(133)
    pcolormesh(reshape(abs.(U - Uexact), mc.n * mc.N, mc.n * mc.N))
    colorbar()
    savefig("poisson_test_$(mpi_size())_$(mc.n).png")
end

if mpi_size()>1
    mpi_finalize()
end


text:
# OpenAI gym environment playing in Julia
# Part of the AGRN.jl package

using AGRN
using PyCall

@pyimport gym

function normalize_obs(obs::Array{Float64})
    (atan.(obs) + (pi/2)) / pi
end

function get_actions(outputs::Array{Float64}, amins::Array{Float64}, amaxs::Array{Float64})
    outputs .* (amaxs .- amins) .+ amins
end

function get_actions(outputs::Array{Float64})
    indmax(outputs)
end

function play_env(grn::GRN, env; render::Bool=false)
    env[:seed](0)
    if render
        env[:render](mode="human")
    end
    ob = env[:reset]()
    discrete_actions = false
    if length(env[:action_space][:shape]) == 0
        discrete_actions = true
    else
        amins = Float64.(env[:action_space][:low])
        amaxs = Float64.(env[:action_space][:high])
    end
    total_reward = 0.0
    done = false
    reward = 0.0

    while ~done
        outputs = step!(grn, normalize_obs(ob))
        if discrete_actions
            action = get_actions(outputs)
        else
            action = get_actions(outputs, amins, amaxs)
        end
        try
            ob, reward, done, _ = env[:step](action)
            if render
                env[:render](mode="human")
            end
        catch
            done = true
        end
        total_reward += reward
    end

    total_reward
end

function get_fitness(env_id::String)
    env = gym.make(env_id)
    nin = length(env[:observation_space][:low])
    discrete_actions = false
    if length(env[:action_space][:shape]) == 0
        nout = env[:action_space][:n]
    else
        nout = length(env[:action_space][:low])
    end
    gym_fit = x->play_env(x, env)
    Fitness(nin, nout, gym_fit)
end

# For running interactively

# Logging.configure(filename="example.log", level=Logging.INFO)
# fitness = get_fitness("MountainCarContinuous-v0")
# config = AGRN.Config()
# maxfit, best = evolve(fitness, config)


text:
# script for generating the patches for a propeller model
using LinearAlgebra
include("../src/SplinePlotting.jl")
plotly()

function getLineProj(R, theta, theta1)
    #gets the coordinates of the points on the tangent to the arc
    pointCenter = zeros(2)
    pointRight = zeros(2)
    pointLeft = zeros(2)

    pointCenter[1] = R*cos(theta)
    pointCenter[2] = R*sin(theta)

    pointRight[1] = (R*sin(theta)+R*cot(theta)*cos(theta))/(tan(theta+theta1)+cot(theta))
    pointRight[2] = tan(theta+theta1)*pointRight[1]

    pointLeft[1] = (R*sin(theta)+R*cot(theta)*cos(theta))/(tan(theta-theta1)+cot(theta))
    pointLeft[2] = tan(theta-theta1)*pointLeft[1]
    return pointLeft, pointRight, pointCenter
end

#dimensions of the hub
radInt = 5 #inner radius
radExt = 6 #outer radius
hubHeight = 10 #height
centerOffset = 0.5 #amount by which the blade interface is lifted from the hub
P_D = 1.4 # P/D pitch distribution

numBlades = 8  #number of blades

bladeThicknessRatio = 0.02 #
bladeHeightRatio = 0.85 #
bladeRotAngle = [-0.8*(pi/2-atan(P_D/pi)),-0.88*(pi/2-atan(P_D/pi)),-0.96*(pi/2-atan(P_D/pi)), -(pi/2-atan(P_D/pi)),-(pi/2-atan(P_D/pi)), -(pi/2-atan(P_D/pi)), -(pi/2-atan(P_D/pi))]
###########BLADE SHAPE DATA#############
coeff1 = 0.2
coeff2 = 3.0
coeff3 = 0.1
coeff4 = 1.2
#top and bottom surface of the blade
coefs_ext = zeros(4,7,3)
coefs_int = zeros(4,7,3)
coefs_ext[:,:,1] = [  0   30   60   90   120  165  170
                      0    20*coeff4        29*coeff4         30*coeff4         24*coeff4        0*coeff4         -45
                      -0.75*coeff1  -0.65*coeff1  -0.50*coeff1  -0.35*coeff1  -0.20*coeff1  0.0*coeff1  0.0*coeff1
                      1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000]

coefs_ext[:,:,2] = [  0   30   60   90   120  165  170
                      -15  -13 -12  -12 -14  -40  -50
                       -0.75*coeff2  -0.65*coeff2  -0.50*coeff2  -0.35*coeff2  -0.20*coeff2  0.0*coeff2  0.0*coeff2
                      1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000]

coefs_ext[:,:,3] = [  0   30   60   90   120  165  170
                      -60  -67 -74  -81 -88  -89  -58
                       -0.75*coeff3  -0.65*coeff3  -0.50*coeff3  -0.35*coeff3  -0.20*coeff3  0.0*coeff3  0.0*coeff3
                      1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000]
#midsurface of the blade
coefs_int[:,:,1] = [ 0   30   60   90   120  165  170
                     0+2  20*coeff4+2  29*coeff4+2  30*coeff4+2  24*coeff4+2  0*coeff4+2  -45+2
                    -0.75*coeff1/2  -0.65*coeff1/2  -0.50*coeff1/2  -0.35*coeff1/2  -0.20*coeff1/2  0.0*coeff1/2  0.0*coeff1/2
                     1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000]

coefs_int[:,:,2] = [ 0   30   60   90   120  165  175
                      -15  -13 -12  -12 -14  -40  -50
                     -0.75*coeff2/2  -0.65*coeff2/2  -0.50*coeff2/2  -0.35*coeff2/2  -0.20*coeff2/2  0.0*coeff2/2  0.0*coeff2/2
                     1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000];

coefs_int[:,:,3] = [ 0   30   60   90   120  165  170
                     -60-2  -67-2 -74-2  -81-2 -88-2  -89-2  -58-2
                     -0.75*coeff3/2  -0.65*coeff3/2  -0.50*coeff3/2  -0.35*coeff3/2  -0.20*coeff3/2  0.0*coeff3/2  0.0*coeff3/2
                     1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000];
#####END BLADE SHAPE DATA######

arcBlade = 2*pi/numBlades
centerBottom  = [0, 0, 0]

bladeHeightBottom = hubHeight*(1-bladeHeightRatio)/2
bladeHeightMiddle = hubHeight*bladeHeightRatio
bladeHeightTop = bladeHeightBottom
bladeBaseArcRatio = 0.5
bladeBaseBottom = bladeHeightBottom*(1-bladeBaseArcRatio)
bladeBaseTop = bladeHeightBottom + bladeHeightMiddle + bladeBaseArcRatio*bladeHeightTop
bladeCenterHeight = hubHeight/2

arcLeft = arcBlade*(1-bladeThicknessRatio)/2
arcRight = arcLeft
arcMiddle = arcBlade*bladeThicknessRatio

#define the arcs
nrbArcLeftExt = nrbcirc(radExt, centerBottom, 0, arcLeft)
nrbArcMiddleExt = nrbcirc(radExt, centerBottom, arcLeft, arcLeft+arcMiddle)
nrbArcRightExt = nrbcirc(radExt, centerBottom, arcLeft+arcMiddle, arcLeft+arcMiddle+arcRight)

nrbArcLeftInt = nrbcirc(radInt, centerBottom, 0, arcLeft)
nrbArcMiddleInt = nrbcirc(radInt, centerBottom, arcLeft, arcLeft+arcMiddle)
nrbArcRightInt = nrbcirc(radInt, centerBottom, arcLeft+arcMiddle, arcLeft+arcMiddle+arcRight)

knotExt = [0,0,0,1/3,1/3,2/3,2/3,1,1,1]
coefsExt = hcat(nrbArcLeftExt.coefs, nrbArcMiddleExt.coefs[:,2:end], nrbArcRightExt.coefs[:,2:end])
nrbArcExt = nrbmak(coefsExt, [knotExt])

knotInt = [0,0,0,1/3,1/3,2/3,2/3,1,1,1]
coefsInt = hcat(nrbArcLeftInt.coefs, nrbArcMiddleInt.coefs[:,2:end], nrbArcRightInt.coefs[:,2:end])
nrbArcInt = nrbmak(coefsInt, [knotInt])

knotThickness = [0,0,1,1]
coefsAnnulus = zeros(4,7,2)
coefsAnnulus[:,:,1] = coefsInt
coefsAnnulus[:,:,2] = coefsExt
nrbAnnulus = nrbmak(coefsAnnulus,[knotInt,knotThickness])

# extrude the annuli along the height
layerSet = [bladeBaseBottom, bladeHeightBottom, bladeCenterHeight, bladeHeightBottom+bladeHeightMiddle,  bladeBaseTop, hubHeight]
knotHeight = [0,0,0, 1/3, 1/3, 2/3, 2/3, 1, 1,1]

# ajust the z-coordinate for each layer
coefsHub = zeros(size(coefsAnnulus)...,length(layerSet)+1...)
coefsHub[:,:,:,1] = coefsAnnulus
for indexSlice = 1:length(layerSet)
    currentSlice = layerSet[indexSlice]
    coefsHub[:,:,:,indexSlice+1] = coefsAnnulus
    coefsHub[3,:,:,indexSlice+1] = currentSlice.*coefsAnnulus[4,:,:]
end
nrbHub = nrbmak(coefsHub,[knotInt, knotThickness, knotHeight])

# % hold on
theta = arcBlade/2
theta1 = theta - arcLeft

pointLeft, pointRight, pointCenter = getLineProj(radExt+centerOffset, theta, theta1)

#define the rotation axis and the rotation matrix
ux = cos(theta)
uy = sin(theta)
uz = 0
rotMatrix = [cos(bladeRotAngle[1])+ux^2*(1-cos(bladeRotAngle[1])) ux*uy*(1-cos(bladeRotAngle[1]))-uz*sin(bladeRotAngle[1]) ux*uz*(1-cos(bladeRotAngle[1]))+uy*sin(bladeRotAngle[1])
    uy*ux*(1-cos(bladeRotAngle[1]))+uz*sin(bladeRotAngle[1]) cos(bladeRotAngle[1])+uy^2*(1-cos(bladeRotAngle[1])) uy*uz*(1-cos(bladeRotAngle[1]))-ux*sin(bladeRotAngle[1])
    uz*ux*(1-cos(bladeRotAngle[1]))-uy*sin(bladeRotAngle[1]) uz*uy*(1-cos(bladeRotAngle[1]))+ux*sin(bladeRotAngle[1]) cos(bladeRotAngle[1])+uz^2*(1-cos(bladeRotAngle[1]))]

bladeInterfacePoints = [pointLeft pointCenter pointRight pointLeft pointCenter pointRight pointLeft pointCenter pointRight
                        ones(1,3)*layerSet[2] ones(1,3)*layerSet[3] ones(1,3)*layerSet[4]]
centerPoint = [pointCenter' layerSet[3]]
bladeInterfacePoints[1,:] = bladeInterfacePoints[1,:] .- centerPoint[1]
bladeInterfacePoints[2,:] = bladeInterfacePoints[2,:] .- centerPoint[2]
bladeInterfacePoints[3,:] = bladeInterfacePoints[3,:] .- centerPoint[3]
bladeInterfacePoints = rotMatrix*bladeInterfacePoints
bladeInterfacePoints[1,:] = bladeInterfacePoints[1,:] .+ centerPoint[1]
bladeInterfacePoints[2,:] = bladeInterfacePoints[2,:] .+ centerPoint[2]
bladeInterfacePoints[3,:] = bladeInterfacePoints[3,:] .+ centerPoint[3]

coefsNewHub = coefsHub

newCoefsLow = [bladeInterfacePoints[:,1:3]; ones(1,3)]
newCoefsMid = [bladeInterfacePoints[:,4:6]; ones(1,3)]
newCoefsTop = [bladeInterfacePoints[:,7:9]; ones(1,3)]
coefsNewHub[:,3:5,2,3] = newCoefsLow
coefsNewHub[:,3:5,2,4] = newCoefsMid
coefsNewHub[:,3:5,2,5] = newCoefsTop

nrbNewHub = nrbmak(coefsNewHub, [knotInt, knotThickness, knotHeight])
# Plots.surface(reuse = false)
# plt = nrbctrlplot(nrbNewHub)
# display(plt)

bladeWidth = norm(pointLeft-pointRight)
bladeHeight = layerSet[4]-layerSet[2]

coefs = coefs_ext
coefs_blade = zeros(size(coefs)...,3...)
scaleFactor = bladeHeight/60
coefs[1:2,:,:] = scaleFactor*coefs[1:2,:,:]
coefs_blade[:,:,:,1] = coefs
coefs[3,:,:] .= bladeWidth
coefs_blade[:,:,:,3] = coefs
coefs = coefs_int
coefs[1:2,:,:] = scaleFactor*coefs[1:2,:,:]
coefs[3,:,:] .= bladeWidth/2
coefs_blade[:,:,:,2] = coefs

numix = 5
numiy = 1
p=2
q=2

# initialize the knot vector
knotU = vcat(zeros(p),collect(LinRange(0,1,numix+1)),ones(p))
knotV = vcat(zeros(q),collect(LinRange(0,1,numiy+1)),ones(q))
knotW = vcat(zeros(q+1), ones(q+1))

# rotate about the x-axis by pi/2
thetaX = pi/2
rotMatrixX = [1  0  0; 0  cos(thetaX)  -sin(thetaX); 0  sin(thetaX)  cos(thetaX)]
for i=1:size(coefs_blade,3)
    for j=1:size(coefs_blade,4)
        coefs_blade[1:3,:,i,j] = rotMatrixX*coefs_blade[1:3,:,i,j]
    end
end

# translate back (in the positive direction of y-axis) by bladeWidth/2
coefs_blade[2,:,:,:] = coefs_blade[2,:,:,:].+bladeWidth/2
# rotate about the z-axis by theta
rotMatrixZ = [cos(theta)  -sin(theta)  0; sin(theta)  cos(theta)  0; 0  0  1]
for i=1:size(coefs_blade,3)
    for j=1:size(coefs_blade,4)
        coefs_blade[1:3,:,i,j] = rotMatrixZ*coefs_blade[1:3,:,i,j]
    end
end

# shift to the center of the target blade location
coefs_blade[3,:,:,:] = coefs_blade[3,:,:,:] .+ bladeHeight/2;
# rotate the blade by bladeRot
for i=1:size(coefs_blade,3)
    for j=1:size(coefs_blade,4)
        for k=1:size(coefs_blade,2)
            global rotMatrix
            R11 = cos(bladeRotAngle[k])+ux^2*(1-cos(bladeRotAngle[k]))
            R12 = ux*uy*(1-cos(bladeRotAngle[k]))-uz*sin(bladeRotAngle[k])
            R13 = ux*uz*(1-cos(bladeRotAngle[k]))+uy*sin(bladeRotAngle[k])
            R21 = uy*ux*(1-cos(bladeRotAngle[k]))+uz*sin(bladeRotAngle[k])
            R22 = cos(bladeRotAngle[k])+uy^2*(1-cos(bladeRotAngle[k]))
            R23 = uy*uz*(1-cos(bladeRotAngle[k]))-ux*sin(bladeRotAngle[k])
            R31 = uz*ux*(1-cos(bladeRotAngle[k]))-uy*sin(bladeRotAngle[k])
            R32 = uz*uy*(1-cos(bladeRotAngle[k]))+ux*sin(bladeRotAngle[k])
            R33 = cos(bladeRotAngle[k])+uz^2*(1-cos(bladeRotAngle[k]))
            rotMatrix = [ R11 R12 R13; R21 R22 R23; R31 R32 R33]
            coefs_blade[1:3,k,i,j] = rotMatrix*coefs_blade[1:3,k,i,j]
        end
    end
end

# shift to the target location
coefs_blade[1,:,:,:] = coefs_blade[1,:,:,:] .+ centerPoint[1]
coefs_blade[2,:,:,:] = coefs_blade[2,:,:,:] .+ centerPoint[2]
coefs_blade[3,:,:,:] = coefs_blade[3,:,:,:] .+ centerPoint[3]
nrbBlade = nrbmak(coefs_blade, [knotU, knotV,knotW]);

# Plots.surface(reuse = false)
# plt = nrbctrlplot(nrbBlade)
# display(plt)

#rotate nrbNewHub and nrbBlade to get the full propeller
Plots.surface(reuse = false)

hubPatches = Array{NURBS, 1}(undef, numBlades)
bladePatches = Array{NURBS, 1}(undef, numBlades)
for i=1:numBlades
    global plt
    global rotMatrixZ
    rotAngle = arcBlade*(i-1);
    rotMatrixZ = [cos(rotAngle) -sin(rotAngle) 0; sin(rotAngle) cos(rotAngle) 0; 0 0 1];
    #rotate and plot the blade
    coefs_blade_current = coefs_blade;
    for j=1:size(coefs_blade,3)
        for k=1:size(coefs_blade,4)
            coefs_blade_current[1:3,:,j,k] = rotMatrixZ*coefs_blade[1:3,:,j,k]
        end
    end
    bladePatches[i] = nrbmak(coefs_blade_current, [knotU, knotV, knotW])
    plt = nrbkntplot(bladePatches[i])
    #plt = nrbplot(bladePatches{i},[20,20,2])

    #rotate and plot the hub
    coefsNewHubCurrent = coefsNewHub;
    for j=1:size(coefsNewHub,3)
        for k=1:size(coefsNewHub,4)
            coefsNewHubCurrent[1:3,:,j,k] = rotMatrixZ*coefsNewHub[1:3,:,j,k]
        end
    end
    hubPatches[i] = nrbmak(coefsNewHubCurrent, [knotInt, knotThickness, knotHeight])
    plt = nrbkntplot(hubPatches[i])
    #plt = nrbplot(hubPatches[i],[20,2,20])
end
display(plt)


text:
Base.@kwdef mutable struct Face
    N = nothing
    NG = nothing
    E = nothing
    F = nothing
    C = nothing
    n = nothing
    a = nothing
    Mid = nothing
    OrientE = nothing
    Dir = nothing
    Ori = nothing
    SF = nothing
    Kite = nothing
    DF = nothing
    DFInv = nothing
    J = nothing
    lon = nothing
    lat = nothing
    P = nothing
end

function Face(Edges,Grid,Pos,Dir,OrientFace=nothing,P=nothing)
  F = Face(;)
  if Edges[1]==0
    return (F,Grid)
  end
  @assert !isnothing(OrientFace)

  F.DF=zeros(0,0);
  nE=size(Edges,2);
  F.F=Pos;
  F.Dir=Dir;
  F.C=zeros(1,2);
  # TODO: check translation
  for iE=1:nE
    Grid.Edges[Edges[iE]].F=[Grid.Edges[Edges[iE]].F Pos];
  end
  #Sort edges
  F.E=zeros(Int,1,nE);
  F.E[1]=Edges[1];
  N2=Grid.Edges[F.E[1]].N[2];
  for iE=2:nE
    for iE1=iE:nE
      if N2==Grid.Edges[Edges[iE1]].N[1]
        F.E[iE]=Edges[iE1];
        N2=Grid.Edges[Edges[iE1]].N[2];
        Edges[iE1]=Edges[iE];
        Edges[iE]=F.E[iE];
        break
      elseif N2==Grid.Edges[Edges[iE1]].N[2]
        F.E[iE]=Edges[iE1];
        N2=Grid.Edges[Edges[iE1]].N[1];
        Edges[iE1]=Edges[iE];
        Edges[iE]=F.E[iE];
        break
      end
    end
  end
  F.a=0;
  F.N=zeros(Int, 1,nE);
  F.N[1,1:2]=Grid.Edges[F.E[1]].N;
  for iE=2:nE-1
    if F.N[1,iE]==Grid.Edges[F.E[iE]].N[1]
      F.N[1,iE+1]=Grid.Edges[F.E[iE]].N[2];
    else
      F.N[1,iE+1]=Grid.Edges[F.E[iE]].N[1];
    end
  end
  if P == nothing
    F.P=zeros(3,size(F.N,2));
    for i=1:size(F.N,2)
      F.P[:,i]=Grid.Nodes[F.N[i]].P';
    end
  else
    F.P=P;
  end
  PT=[0, 0, 0];
  for i=1:nE-1
    #P=P+cross(Grid.Nodes(F.N(i)).P,Grid.Nodes(F.N(i+1)).P);
    PT=PT+cross(F.P[:,i],F.P[:,i+1]);
  end
  #P=P+cross(Grid.Nodes(F.N(nE)).P,Grid.Nodes(F.N[1]).P);
  PT=PT+cross(F.P[:,nE],F.P[:,1]);
  F.a=0.5*norm(PT);
  F.Mid=zeros(3);
  for i=1:nE
    #F.Mid=F.Mid+Grid.Nodes(F.N(i)).P;
    F.Mid=F.Mid+F.P[:,i];
  end
  F.Mid=F.Mid/nE;

  NumE=size(Edges,2);
  #F.n=cross(Grid.Nodes(F.N(NumE)).P,Grid.Nodes(F.N[1]).P);
  F.n=cross(F.P[:,NumE],F.P[:,1]);
  for i=1:NumE-1
    #F.n=F.n+cross(Grid.Nodes(F.N(i)).P,Grid.Nodes(F.N(i+1)).P);
    F.n=F.n+cross(F.P[:,i],F.P[:,i+1]);
  end
  F.n=F.n/norm(F.n);
  if OrientFace(F.n',F.Mid)<0 # TODO: check translation with Oswald
    #Change Orientation
    NTemp=copy(F.N);
    ETemp=copy(F.E);
    PTemp=copy(F.P);
    for i=1:nE
      F.N[i]=NTemp[nE-i+1];
      F.P[:,i]=PTemp[:,nE-i+1];
    end
    for i=1:nE-1
      F.E[i]=ETemp[nE-i];
    end
    F.n=-F.n;
  end
  return (F,Grid)
end