diff --git a/.gitignore b/.gitignore index 7269e4bc..7cb5316c 100644 --- a/.gitignore +++ b/.gitignore @@ -39,3 +39,7 @@ _deps #vs code settings .vscode .DS_Store + +# fism files +srcPython/FISM*.nc +srcPython/fism*.txt \ No newline at end of file diff --git a/doc/internals/fism.md b/doc/internals/fism.md new file mode 100644 index 00000000..5b1322e8 --- /dev/null +++ b/doc/internals/fism.md @@ -0,0 +1,42 @@ +# FISM-2 + +FISM-2 can be used as a EUV model. + +Needs the FISM-2 files automatically made by `srcPython/fism.py`. + +FISM-2 contains the binned flux already, however Aether must still be provided with an +EUV csv file containing the cross-sections. + +The different fism models available in fism.py contain different numbers of bins, so +the euv file must be different. + +| Model | number of bins | euv file | +| :--- | :-------------: | -------: | +| HFG | 23 | euv_solomon.csv | +| Solomon | 23 | euv_solomon.csv | +| NEUVAC | 37/59 | euv.csv / euv_59.csv | +| EUVAC | 37 | euv.csv | + +The input format, when using fism data: + + "Euv" : { + "doUse" : true, + "Model" : "fism", + "File" : "UA/inputs/euv_59.csv", + "fismFile": "fism2_file_59.txt", + "IncludePhotoElectrons" : true, + "HeatingEfficiency" : 0.05, + "dt" : 60.0 + }, + +To generate FISM-2 irradiances between two dates, one may call fism.py like so: + +srcPython/fism.py 20110319 20110321 -b neuvac + +Note that the optional argument '-b' defaults to the binning scheme of the NEUVAC model, +which employs 59 bins. To use the 37 bins of the EUVAC model, the argument should be +'euvac', while to use the 23 bins used by the HFG model, the argument should be 'solomon'. + +fism.py should always be run before Aether is run using the FISM-2 model. Even though the +entire FISM-2 irradiances are stored in the repository, Aether will need a separate .csv +file containing the temporal subset of FISM-2 irradiances in the desired binning scheme. diff --git a/include/euv.h b/include/euv.h index a9dd826d..dca2825a 100644 --- a/include/euv.h +++ b/include/euv.h @@ -94,6 +94,11 @@ class Euv { /// NEUVAC model intercept: std::vector neuvac_int; + // To avoid having to start from 0 each iteration: + int fism_prev_index = 0; + // Declare this so it is not passed between function: + index_file_output_struct fismData; + // -------------------------------------------------------------------- // Functions: @@ -109,13 +114,26 @@ class Euv { **/ bool euvac(Times time, Indices indices); + /********************************************************************** + \brief Compute the EUV spectrum given F107 and F107a + \param time The times within the model (dt is needed) + \param indices Need the F107 and F107a + **/ + bool solomon_hfg(Times time, Indices indices); + /********************************************************************** - \brief Compute the EUV spectrum given F107 and F107a - \param time The times within the model (dt is needed) - \param indices Need the F107 and F107a - **/ - bool solomon_hfg(Times time, Indices indices); + \brief returns the FISM spectrum for a given time + + Unlike the other EUV models ([N]EUVAC, Solomon, etc.), the spectrum + is read from a file (stored in fismData). This does the same thing + as get_index, however FISM is not stored in Indices since it can + have variable # of bins + + \param time The times within the model (dt is needed) + **/ + bool get_fism(Times time); + /********************************************************************** \brief Compute the EUV spectrum given F107 and F107a (new version) \param time The times within the model (dt is needed) @@ -161,6 +179,17 @@ class Euv { **/ bool read_file(); + /********************************************************************** + \brief Read in the FISM file + + Read in the CSV file with FISM data. This can be made with + srcPython/fism.py. The data are read into a index_file_output_struct, + where each row is one time, and each col is a "variable". These should + match the number of bins in the provided EUV file. + **/ + index_file_output_struct read_fism(std::string fism_filename); + + /********************************************************************** \brief Interprets the EUV CSV rows and returns the relevant row diff --git a/include/inputs.h b/include/inputs.h index 7c48dbf1..793859ed 100644 --- a/include/inputs.h +++ b/include/inputs.h @@ -228,6 +228,12 @@ class Inputs { \brief returns settings[" \param **/ + std::string get_euv_fismfile(); + + /********************************************************************** + \brief returns settings[" + \param + **/ bool get_euv_douse(); /********************************************************************** diff --git a/share/run/UA/inputs/euv_59_v2.csv b/share/run/UA/inputs/euv_59_v2.csv new file mode 100644 index 00000000..7f807952 --- /dev/null +++ b/share/run/UA/inputs/euv_59_v2.csv @@ -0,0 +1,43 @@ + Short,,wave,1.0000E+00,Angstroms, 1.0000E+00, 2.0000E+00, 4.0000E+00, 8.0000E+00, 1.6000E+01, 2.3000E+01, 3.2000E+01, 5.0000E+01, 1.0000E+02, 1.5000E+02, 2.0000E+02, 2.5630E+02, 2.8415E+02, 2.5000E+02, 3.0331E+02, 3.0378E+02, 3.0000E+02, 3.6807E+02, 3.5000E+02, 4.0000E+02, 4.6522E+02, 4.5000E+02, 5.0000E+02, 5.5437E+02, 5.8433E+02, 5.5000E+02, 6.0976E+02, 6.2973E+02, 6.0000E+02, 6.5000E+02, 7.0331E+02, 7.0000E+02, 7.6515E+02, 7.7041E+02, 7.8936E+02, 7.5000E+02, 8.0000E+02, 8.5000E+02, 9.0000E+02, 9.7702E+02, 9.5000E+02, 1.0257E+03, 1.0319E+03, 1.0000E+03, 1.0500E+03, 1.1000E+03, 1.1500E+03, 1.2157E+03, 1.2000E+03, 1.2500E+03, 1.3000E+03, 1.3500E+03, 1.4000E+03, 1.4500E+03, 1.5000E+03, 1.5500E+03, 1.6000E+03, 1.6500E+03, 1.7000E+03, from GITM + Long,,wave,1.0000E+00,Angstroms, 2.0000E+00, 4.0000E+00, 8.0000E+00, 1.6000E+01, 2.3000E+01, 3.2000E+01, 5.0000E+01, 1.0000E+02, 1.5000E+02, 2.0000E+02, 2.5000E+02, 2.5630E+02, 2.8415E+02, 3.0000E+02, 3.0331E+02, 3.0378E+02, 3.5000E+02, 3.6807E+02, 4.0000E+02, 4.5000E+02, 4.6522E+02, 5.0000E+02, 5.5000E+02, 5.5437E+02, 5.8433E+02, 6.0000E+02, 6.0976E+02, 6.2973E+02, 6.5000E+02, 7.0000E+02, 7.0331E+02, 7.5000E+02, 7.6515E+02, 7.7041E+02, 7.8936E+02, 8.0000E+02, 8.5000E+02, 9.0000E+02, 9.5000E+02, 9.7702E+02, 1.0000E+03, 1.0257E+03, 1.0319E+03, 1.0500E+03, 1.1000E+03, 1.1500E+03, 1.2000E+03, 1.2157E+03, 1.2500E+03, 1.3000E+03, 1.3500E+03, 1.4000E+03, 1.4500E+03, 1.5000E+03, 1.5500E+03, 1.6000E+03, 1.6500E+03, 1.7000E+03, 1.7500E+03, from GITM + F74113,,,1.0000E+09,/cm2/s, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 1.2000E+00, 4.5000E-01, 4.8000E+00, 3.1000E+00, 4.6000E-01, 2.1000E-01, 1.6790E+00, 8.0000E-01, 6.9000E+00, 9.6500E-01, 6.5000E-01, 3.1400E-01, 3.8300E-01, 2.9000E-01, 2.8500E-01, 4.5200E-01, 7.2000E-01, 1.2700E+00, 3.5700E-01, 5.3000E-01, 1.5900E+00, 3.4200E-01, 2.3000E-01, 3.6000E-01, 1.4100E-01, 1.7000E-01, 2.6000E-01, 7.0200E-01, 7.5800E-01, 1.6250E+00, 3.5370E+00, 3.0000E+00, 4.4000E+00, 1.4750E+00, 3.5000E+00, 2.1000E+00, 2.4670E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + AFAC,,,1.0000E+00,, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 1.0017E-02, 7.1250E-03, 1.3375E-02, 1.9450E-02, 2.7750E-03, 1.3768E-01, 2.6467E-02, 2.5000E-02, 3.3333E-03, 2.2450E-02, 6.5917E-03, 3.6542E-02, 7.4083E-03, 7.4917E-03, 2.0225E-02, 8.7583E-03, 3.2667E-03, 5.1583E-03, 3.6583E-03, 1.6175E-02, 3.3250E-03, 1.1800E-02, 4.2667E-03, 3.0417E-03, 4.7500E-03, 3.8500E-03, 1.2808E-02, 3.2750E-03, 4.7667E-03, 4.8167E-03, 5.6750E-03, 4.9833E-03, 3.9417E-03, 4.4167E-03, 5.1833E-03, 5.2833E-03, 4.3750E-03, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + NEUV_S1,,,1.0000e+00,slope, 4.7721e-01, 4.8836e-01, 5.1112e-01, 4.6218e-01, 4.4626e-01, 4.5621e-01, 4.6444e-01, 2.5497e-01, 4.5348e-01, 2.7143e-01, 1.0616e+00, 4.3086e-01, 4.8995e-01, 4.3046e-01, 1.9877e-01, 4.5349e-01, 5.1979e-01, 4.5590e-01, 4.0438e-01, 4.4541e-01, 4.5721e-01, 5.2010e-01, 4.5938e-01, 4.5614e-01, 4.9421e-01, 4.6293e-01, 4.6222e-01, 4.5820e-01, 4.4724e-01, 4.5370e-01, 4.5756e-01, 4.6456e-01, 4.5873e-01, 4.5651e-01, 4.5443e-01, 4.4956e-01, 4.1966e-01, 5.3686e-01, 4.3159e-01, 3.6366e-01, 4.7233e-01, 4.7194e-01, 4.5386e-01, 4.5133e-01, 4.6587e-01, 4.4840e-01, 5.3661e-01, 1.0786e-01, 2.6182e-01, 4.4807e-01, 4.2671e-01, 4.9602e-01, 4.4946e-01, 5.0319e-01, 4.3481e-01, 4.9051e-01, 1.5122e-01, 3.8279e-01, 7.0315e-01,from GITM + NEUV_S2,,,1.0000e+00,slope, -4.7721e-01, -4.8836e-01, -5.1112e-01, -4.6218e-01, -4.4626e-01, -4.5620e-01, -4.6443e-01, -2.5496e-01, -4.5348e-01, -2.7142e-01, -1.0616e+00, -4.3086e-01, -4.8995e-01, -4.3045e-01, -1.9877e-01, -4.5349e-01, -5.1978e-01, -4.5590e-01, -4.0438e-01, -4.4540e-01, -4.5721e-01, -5.2010e-01, -4.5938e-01, -4.5614e-01, -4.9421e-01, -4.6293e-01, -4.6222e-01, -4.5820e-01, -4.4724e-01, -4.5370e-01, -4.5756e-01, -4.6456e-01, -4.5873e-01, -4.5651e-01, -4.5443e-01, -4.4956e-01, -4.1966e-01, -5.3685e-01, -4.3158e-01, -3.6366e-01, -4.7233e-01, -4.7194e-01, -4.5386e-01, -4.5133e-01, -4.6587e-01, -4.4840e-01, -5.3661e-01, -1.0781e-01, -2.6182e-01, -4.4807e-01, -4.2671e-01, -4.9602e-01, -4.4946e-01, -5.0318e-01, -4.3481e-01, -4.9051e-01, -1.5121e-01, -3.8278e-01, -7.0314e-01,from GITM + NEUV_S3,,,1.0000e+00,slope, 4.7721e-01, 4.8836e-01, 5.1112e-01, 4.6218e-01, 4.4626e-01, 4.5621e-01, 4.6444e-01, 2.5496e-01, 4.5348e-01, 2.7143e-01, 1.0616e+00, 4.3086e-01, 4.8995e-01, 4.3045e-01, 1.9877e-01, 4.5349e-01, 5.1979e-01, 4.5590e-01, 4.0438e-01, 4.4541e-01, 4.5721e-01, 5.2010e-01, 4.5938e-01, 4.5614e-01, 4.9421e-01, 4.6293e-01, 4.6222e-01, 4.5820e-01, 4.4724e-01, 4.5370e-01, 4.5756e-01, 4.6456e-01, 4.5873e-01, 4.5651e-01, 4.5443e-01, 4.4956e-01, 4.1966e-01, 5.3685e-01, 4.3158e-01, 3.6366e-01, 4.7233e-01, 4.7194e-01, 4.5386e-01, 4.5133e-01, 4.6587e-01, 4.4840e-01, 5.3661e-01, 1.0783e-01, 2.6182e-01, 4.4807e-01, 4.2671e-01, 4.9602e-01, 4.4946e-01, 5.0319e-01, 4.3481e-01, 4.9051e-01, 1.5121e-01, 3.8278e-01, 7.0314e-01,from GITM + NEUV_l1,,,1.0000e+00,ints, 6.3336e-12, -1.7586e-09, -5.8173e-07, -2.7526e-05, -7.9475e-05, -2.5214e-05, -7.7868e-05, -1.8471e-04, -1.0315e-05, -4.5460e-06, -1.7387e-04, -2.8295e-06, -9.3578e-05, -8.0510e-05, 2.5654e-04, 9.7116e-07, -4.2748e-05, 1.8657e-05, -4.2105e-05, 1.1967e-05, 1.3408e-05, -3.2775e-06, 7.3782e-06, 2.1466e-05, 2.6201e-05, 8.4018e-06, -1.1750e-06, 4.4671e-05, 2.3090e-07, 1.1496e-05, 1.0248e-05, 1.1219e-05, 1.0442e-05, 1.1598e-05, 1.0931e-05, 3.5491e-05, 5.7380e-05, 7.2432e-05, 6.6057e-05, 9.0619e-05, 3.0050e-05, 3.3145e-05, 2.8926e-05, 3.9467e-05, 3.6713e-05, 4.8318e-05, 1.5926e-04, 3.6572e-03, 3.2218e-04, 8.4120e-05, 2.6237e-04, 1.4758e-04, 1.9098e-04, 2.9104e-04, 4.6908e-04, 7.1953e-04, 1.0507e-03, 1.9508e-03, 3.2611e-03,from GITM + NEUV_P1,,,1.0000e+00,powers, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 9.9997e-01, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00,from GITM + NEUV_P2,,,1.0000e+00,powers, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00,from GITM + O2,O2,abs,1.0000E-22,m2, 4.0000E-04, 4.0000E-03, 2.4000E-02, 1.4000E-01, 1.0200E+00, 1.0000E-01, 3.2000E-01, 1.1800E+00, 4.0000E+00, 7.1000E+00, 1.0600E+01, 1.3200E+01, 1.5700E+01, 1.5100E+01, 1.6800E+01, 1.6800E+01, 1.7190E+01, 1.8400E+01, 1.8170E+01, 1.9390E+01, 2.0400E+01, 2.1590E+01, 2.4060E+01, 2.5590E+01, 2.2000E+01, 2.5040E+01, 2.6100E+01, 2.5800E+01, 2.6020E+01, 2.6270E+01, 2.5000E+01, 2.9050E+01, 2.1960E+01, 2.5180E+01, 2.6660E+01, 2.7090E+01, 2.0870E+01, 9.8500E+00, 1.5540E+01, 4.0000E+00, 1.6530E+01, 1.6000E+00, 1.0000E+00, 1.1000E+00, 1.0000E+00, 1.0000E-01, 3.0000E-01, 1.0000E-02, 3.0000E+00, 3.0000E-01, 2.2000E+00, 1.2000E+01, 1.5000E+01, 1.3000E+01, 1.0000E+01, 6.0000E+00, 3.4000E+00, 1.5000E+00, 5.0000E-01, from GITM + O,O,abs,1.0000E-22,m2, 2.0000E-04, 2.0000E-03, 1.2000E-02, 7.0000E-02, 5.1000E-01, 5.0000E-02, 1.6000E-01, 5.9000E-01, 1.6000E+00, 2.9000E+00, 5.3000E+00, 6.0500E+00, 7.1300E+00, 6.6100E+00, 7.6800E+00, 7.7000E+00, 8.6700E+00, 9.9500E+00, 9.6400E+00, 1.1210E+01, 1.1250E+01, 1.1640E+01, 1.1910E+01, 1.2130E+01, 1.2170E+01, 1.1900E+01, 1.2230E+01, 1.2220E+01, 1.2210E+01, 1.0040E+01, 1.1350E+01, 8.0000E+00, 4.1800E+00, 4.1800E+00, 4.2800E+00, 4.2300E+00, 4.3800E+00, 4.1800E+00, 2.1200E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + N2,N2,abs,1.0000E-22,m2, 3.0000E-04, 3.0000E-03, 1.5000E-02, 9.0000E-02, 4.8000E-01, 1.1600E+00, 2.4000E-01, 6.0000E-01, 1.9000E+00, 4.4000E+00, 8.0000E+00, 9.7000E+00, 1.0600E+01, 1.0300E+01, 1.1600E+01, 1.1600E+01, 1.3000E+01, 1.8000E+01, 1.7510E+01, 2.1070E+01, 2.1800E+01, 2.1850E+01, 2.4530E+01, 2.4690E+01, 2.3200E+01, 2.2380E+01, 2.3100E+01, 2.3200E+01, 2.3220E+01, 2.9750E+01, 2.6300E+01, 3.0940E+01, 3.5460E+01, 2.6880E+01, 1.9260E+01, 3.0710E+01, 1.5050E+01, 4.6630E+01, 1.6990E+01, 7.0000E-01, 3.6160E+01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + CO2,CO2,abs,1.0000E-22,m2, 3.5800E-03, 3.5800E-03, 3.5800E-03, 3.5800E-01, 3.5800E-01, 3.5800E-01, 3.5800E-01, 1.5500E+00, 4.6200E+00, 9.0900E+00, 1.4360E+01, 1.6510E+01, 1.9020E+01, 1.7520E+01, 2.1490E+01, 2.1590E+01, 2.3570E+01, 2.5270E+01, 2.4870E+01, 2.8270E+01, 2.9530E+01, 3.0250E+01, 3.1490E+01, 3.3200E+01, 3.4200E+01, 3.4910E+01, 3.5300E+01, 3.4300E+01, 3.4450E+01, 3.3700E+01, 2.3520E+01, 3.2830E+01, 9.3840E+01, 6.1940E+01, 2.6490E+01, 3.9830E+01, 1.3980E+01, 4.4670E+01, 5.2080E+01, 4.2870E+01, 5.0310E+01, 1.5100E+01, 1.4200E+01, 1.8240E+01, 1.7400E+01, 4.0800E+01, 8.8200E-01, 4.9600E-02, 8.1000E-02, 3.7300E-01, 7.3900E-01, 6.0700E-01, 5.2400E-01, 5.4400E-01, 4.3100E-01, 2.5800E-01, 1.2600E-01, 4.8000E-02, 1.6000E-02, from GITM + CO,CO,abs,1.0000E-22,m2, 4.1700E-03, 4.1700E-03, 4.1700E-03, 4.1700E-01, 4.1700E-01, 4.1700E-01, 4.1700E-01, 8.7000E-01, 2.3900E+00, 4.6700E+00, 7.0100E+00, 8.6100E+00, 1.0540E+01, 9.4200E+00, 1.1870E+01, 1.1900E+01, 1.3440E+01, 1.5260E+01, 1.4960E+01, 1.7960E+01, 2.0170E+01, 2.0570E+01, 2.1090E+01, 2.1620E+01, 2.2000E+01, 2.1910E+01, 2.2100E+01, 2.2030E+01, 2.1920E+01, 2.1040E+01, 2.3850E+01, 2.5200E+01, 2.6280E+01, 1.5260E+01, 3.3130E+01, 2.0540E+01, 2.2610E+01, 3.6980E+01, 5.0320E+01, 2.8500E+01, 5.2830E+01, 1.3900E+00, 1.3900E+00, 8.5700E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + CH4,CH4,abs,1.0000E-22,m2, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 1.6430E-01, 4.2100E-01, 4.5350E-01, 2.0400E-01, 5.9300E-01, 1.4960E+00, 2.7940E+00, 3.8570E+00, 5.0530E+00, 4.3600E+00, 6.0330E+00, 6.0590E+00, 7.8290E+00, 1.0165E+01, 9.7760E+00, 1.4701E+01, 1.8770E+01, 2.1449E+01, 2.4644E+01, 2.7924E+01, 3.1052E+01, 3.0697E+01, 3.3178E+01, 3.5276E+01, 3.4990E+01, 3.9280E+01, 4.1069E+01, 4.2927E+01, 4.5458E+01, 4.5716E+01, 4.6472E+01, 4.5921E+01, 4.8327E+01, 4.8968E+01, 4.8001E+01, 4.1154E+01, 3.8192E+01, 3.2700E+01, 3.0121E+01, 2.9108E+01, 2.8400E+01, 1.8000E+01, 1.9200E+01, 1.7860E+01, 1.8318E+01, 1.9068E+01, 1.2826E+01, 3.2898E+00, 1.2600E-01, 7.9900E-04, 1.4000E-05, 7.0000E-06, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + H2,H2,abs,1.0000E-22,m2, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 1.0000E-02, 8.0000E-02, 2.1000E-01, 4.3000E-01, 6.0000E-01, 8.4000E-01, 7.3000E-01, 1.0200E+00, 1.0200E+00, 1.4200E+00, 1.9400E+00, 1.9000E+00, 3.0300E+00, 3.8700E+00, 4.5000E+00, 5.3600E+00, 6.1700E+00, 7.0200E+00, 6.8600E+00, 7.8100E+00, 8.4600E+00, 8.4500E+00, 9.9000E+00, 1.0730E+01, 1.1370E+01, 1.0760E+01, 8.6400E+00, 7.3400E+00, 8.7500E+00, 8.2500E+00, 4.8000E-01, 1.9000E-01, 0.0000E+00, 5.0000E-02, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + HCN,HCN,abs,1.0000E-22,m2, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + He,He,abs,1.0000E-22,m2, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 1.4400E-01, 4.7900E-01, 1.1570E+00, 1.6010E+00, 2.1210E+00, 2.5950E+00, 2.3200E+00, 2.9530E+00, 2.9620E+00, 3.5440E+00, 4.2680E+00, 4.1420E+00, 5.4470E+00, 6.5630E+00, 7.2080E+00, 9.5800E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + NO,NO,abs,1.0000E-22,m2, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 2.4000E+01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 2.4000E+01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 1.0000E+01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 2.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + O2,O2+,ion,1.0000E-22,m2, 4.0000E-04, 4.0000E-03, 2.4000E-02, 1.4000E-01, 1.0200E+00, 1.0000E-01, 3.2000E-01, 1.1800E+00, 4.0000E+00, 7.1000E+00, 1.0600E+01, 1.3200E+01, 1.5700E+01, 1.5100E+01, 1.6800E+01, 1.6800E+01, 1.7190E+01, 1.8400E+01, 1.8170E+01, 1.9390E+01, 2.0400E+01, 2.1590E+01, 2.4060E+01, 2.5590E+01, 2.2000E+01, 2.5040E+01, 2.6100E+01, 2.5800E+01, 2.5940E+01, 2.2050E+01, 2.3000E+01, 2.3810E+01, 8.5900E+00, 9.6900E+00, 1.1050E+01, 9.3900E+00, 6.1200E+00, 4.6900E+00, 9.3400E+00, 2.5000E+00, 1.2220E+01, 1.0000E+00, 0.0000E+00, 2.7000E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + O,O+,ion,1.0000E-22,m2, 6.0000E-05, 6.0000E-04, 4.0000E-03, 2.0000E-02, 1.5000E-02, 1.5000E-02, 5.0000E-02, 1.8000E-01, 4.6000E-01, 7.8000E-01, 1.3800E+00, 1.5100E+00, 1.7800E+00, 1.6500E+00, 1.9200E+00, 1.9300E+00, 2.2500E+00, 2.5900E+00, 2.5100E+00, 3.0300E+00, 3.1500E+00, 3.2600E+00, 3.4500E+00, 3.5200E+00, 3.5300E+00, 3.4500E+00, 3.6700E+00, 3.7900E+00, 3.7800E+00, 4.0100E+00, 4.9100E+00, 4.2000E+00, 4.1800E+00, 4.1800E+00, 4.2800E+00, 4.2300E+00, 4.3800E+00, 4.1800E+00, 2.1200E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + O,O+2D,ion,1.0000E-22,m2, 6.0000E-05, 6.0000E-04, 4.0000E-03, 2.0000E-02, 1.5000E-02, 1.5000E-02, 5.0000E-02, 1.9000E-01, 5.1000E-01, 9.9000E-01, 1.8600E+00, 2.1200E+00, 2.6400E+00, 2.3800E+00, 2.8400E+00, 2.8500E+00, 3.4700E+00, 3.9800E+00, 3.8600E+00, 4.7100E+00, 5.0600E+00, 5.2400E+00, 5.3600E+00, 5.4600E+00, 5.4800E+00, 5.3600E+00, 5.5000E+00, 5.5000E+00, 5.4900E+00, 5.5200E+00, 6.4400E+00, 3.8000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + O,O+2P,ion,1.0000E-22,m2, 7.0000E-05, 7.0000E-04, 4.0000E-03, 3.0000E-02, 2.0000E-02, 2.0000E-02, 6.0000E-02, 2.2000E-01, 6.2000E-01, 1.1300E+00, 2.0700E+00, 2.4200E+00, 2.7100E+00, 2.5800E+00, 2.9200E+00, 2.9300E+00, 2.9500E+00, 3.3800E+00, 3.2800E+00, 3.4800E+00, 3.0400E+00, 3.1400E+00, 3.1000E+00, 3.1500E+00, 3.1600E+00, 3.0900E+00, 3.0600E+00, 2.9300E+00, 2.9300E+00, 5.0000E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + N2,N2+,ion,1.0000E-22,m2, 3.0000E-04, 3.0000E-03, 1.5000E-02, 9.0000E-02, 4.8000E-01, 1.1600E+00, 2.4000E-01, 6.0000E-01, 1.9000E+00, 4.4000E+00, 8.0000E+00, 9.7000E+00, 1.0600E+01, 1.0300E+01, 1.1600E+01, 1.1600E+01, 1.3000E+01, 1.8000E+01, 1.7510E+01, 2.1070E+01, 2.1800E+01, 2.1850E+01, 2.4530E+01, 2.4690E+01, 2.3200E+01, 2.2380E+01, 2.3100E+01, 2.3200E+01, 2.3220E+01, 2.5060E+01, 2.3000E+01, 2.3200E+01, 2.3770E+01, 1.8390E+01, 1.0180E+01, 1.6750E+01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + N,N+,ion,1.0000E-22,m2, 0.0000E+00, 0.0000E+00, 1.0000E-01, 5.0000E-01, 1.0000E+00, 1.0000E+00, 1.0000E+00, 2.0000E+00, 2.5000E+00, 3.5000E+00, 4.0000E+00, 5.0000E+00, 5.0000E+00, 6.0000E+00, 6.0000E+00, 6.5000E+00, 8.0000E+00, 7.0000E+00, 1.0000E+01, 1.0000E+01, 1.0000E+01, 1.1000E+01, 1.1500E+01, 1.2000E+01, 1.1000E+01, 1.2000E+01, 1.2000E+01, 1.2000E+01, 1.2000E+01, 1.2000E+01, 1.1000E+01, 1.1000E+01, 1.1000E+01, 1.0000E+01, 1.0000E+01, 1.0000E+01, 1.0000E+01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + CO2,CO2+,ion,1.0000E-22,m2, 1.5500E+00, 1.5500E+00, 1.5500E+00, 1.5500E+00, 1.5500E+00, 1.5500E+00, 1.5500E+00, 1.5500E+00, 4.6200E+00, 9.0900E+00, 1.4320E+01, 1.6110E+01, 1.8600E+01, 1.7140E+01, 2.1390E+01, 2.1440E+01, 2.3630E+01, 2.5560E+01, 2.5520E+01, 2.7170E+01, 2.8760E+01, 3.0680E+01, 3.2602E+01, 3.3210E+01, 3.3860E+01, 3.4960E+01, 3.5300E+01, 3.4300E+01, 3.4570E+01, 3.2290E+01, 2.0860E+01, 2.7490E+01, 8.6320E+01, 5.1770E+01, 2.1680E+01, 3.4090E+01, 1.0930E+01, 7.1400E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + CO,CO+,ion,1.0000E-22,m2, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + CH4,CH4+,ion,1.0000E-22,m2, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 2.0000E-01, 5.9000E-01, 1.5000E+00, 2.7900E+00, 3.8600E+00, 5.0500E+00, 4.3600E+00, 6.0300E+00, 6.0600E+00, 7.8300E+00, 1.0170E+01, 9.7800E+00, 1.4700E+01, 1.8770E+01, 2.1450E+01, 2.4640E+01, 2.7920E+01, 3.1050E+01, 3.0700E+01, 3.3180E+01, 3.5280E+01, 3.4990E+01, 3.9280E+01, 4.1070E+01, 4.2930E+01, 4.4800E+01, 4.4800E+01, 4.4610E+01, 4.4690E+01, 4.0280E+01, 2.5530E+01, 1.3860E+01, 1.4000E-01, 4.8000E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + H2,H2+,ion,1.0000E-22,m2, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 1.0000E-02, 8.0000E-02, 2.0000E-01, 4.0000E-01, 5.5000E-01, 7.5000E-01, 6.5000E-01, 9.0000E-01, 9.0000E-01, 1.3000E+00, 1.7800E+00, 1.7400E+00, 2.8900E+00, 3.7800E+00, 4.0500E+00, 5.2500E+00, 6.0500E+00, 6.9000E+00, 6.7400E+00, 7.6700E+00, 8.3000E+00, 8.2900E+00, 9.7000E+00, 1.0730E+01, 9.7600E+00, 8.6200E+00, 7.0700E+00, 5.0700E+00, 6.6300E+00, 9.0000E-02, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + HCN,HCN+,ion,1.0000E-22,m2, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + NO,NO+,ion,1.0000E-22,m2, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 2.4000E+01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 2.4000E+01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 1.0000E+01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 2.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + O2,O2+,pei,1.0000E+00,, 1.3469E+02, 1.3469E+02, 3.2212E+01, 1.3309E+01, 3.9615E+01, 3.9615E+01, 2.8340E+00, 1.0920E+00, 1.0920E+00, 4.1600E-01, 1.8400E-01, 1.8400E-01, 1.8400E-01, 1.8400E-01, 9.0000E-02, 9.0000E-02, 2.4000E-02, 2.4000E-02, 2.4000E-02, 2.4000E-02, 2.4000E-02, 2.4000E-02, 2.4000E-02, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + O2,O+,pei,1.0000E+00,, 7.6136E+01, 7.6136E+01, 1.7944E+01, 6.9810E+00, 2.0338E+01, 2.0338E+01, 1.4370E+00, 5.2100E-01, 5.2100E-01, 1.6300E-01, 5.2000E-02, 5.2000E-02, 5.2000E-02, 5.2000E-02, 1.4000E-02, 1.4000E-02, 1.0000E-03, 1.0000E-03, 1.0000E-03, 1.0000E-03, 1.0000E-03, 1.0000E-03, 1.0000E-03, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + N2,N2+,pei,1.0000E+00,, 2.6399E+02, 2.6399E+02, 6.2570E+01, 2.5213E+01, 8.5400E+00, 8.5400E+00, 6.1420E+00, 2.2880E+00, 2.2880E+00, 7.8600E-01, 3.2400E-01, 3.2400E-01, 3.2400E-01, 3.2400E-01, 1.6900E-01, 1.6900E-01, 3.1000E-02, 3.1000E-02, 3.1000E-02, 3.1000E-02, 3.1000E-02, 3.1000E-02, 3.1000E-02, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + N2,N+,pei,1.0000E+00,, 7.8674E+01, 7.8674E+01, 1.8310E+01, 6.9480E+00, 2.2950E+00, 2.2950E+00, 1.6470E+00, 1.6470E+00, 5.7100E-01, 1.4600E-01, 3.7000E-02, 3.7000E-02, 3.7000E-02, 3.7000E-02, 8.0000E-03, 8.0000E-03, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + CH4,CH4+,pei,1.0000E+00,, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + CH4,CH3+,pei,1.0000E+00,, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + O,O+,pei,1.0000E+00,, 8.1240E+01, 8.1240E+01, 1.8896E+01, 9.4250E+00, 2.8622E+01, 2.8622E+01, 2.0190E+00, 9.0200E-01, 9.0200E-01, 4.7000E-01, 3.2500E-01, 3.2500E-01, 3.2500E-01, 3.2500E-01, 2.0900E-01, 2.0900E-01, 8.4000E-02, 8.4000E-02, 8.4000E-02, 8.4000E-02, 8.4000E-02, 8.4000E-02, 8.4000E-02, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + O,O+2D,pei,1.0000E+00,, 8.8526E+01, 8.8526E+01, 2.0691E+01, 9.3650E+00, 2.8199E+01, 2.8199E+01, 1.9620E+00, 8.5300E-01, 8.5300E-01, 4.1800E-01, 2.5300E-01, 2.5300E-01, 2.5300E-01, 2.5300E-01, 1.4800E-01, 1.4800E-01, 3.4000E-02, 3.4000E-02, 3.4000E-02, 3.4000E-02, 3.4000E-02, 3.4000E-02, 3.4000E-02, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + O,O+2P,pei,1.0000E+00,, 4.7358E+01, 4.7358E+01, 1.1007E+01, 4.7720E+00, 1.4556E+01, 1.4556E+01, 1.0140E+00, 4.3600E-01, 4.3600E-01, 2.0300E-01, 1.1600E-01, 1.1600E-01, 1.1600E-01, 1.1600E-01, 6.1000E-02, 6.1000E-02, 9.0000E-03, 9.0000E-03, 9.0000E-03, 9.0000E-03, 9.0000E-03, 9.0000E-03, 9.0000E-03, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + N2,N_4S,ped,1.0000E+00,, 2.4500E+02, 2.4500E+02, 5.2052E+01, 2.5255E+01, 9.0490E+00, 9.0490E+00, 6.5320E+00, 2.9090E+00, 2.9090E+00, 1.3710E+00, 7.6400E-01, 7.6400E-01, 7.6400E-01, 7.6400E-01, 5.1500E-01, 5.1500E-01, 1.5700E-01, 1.5700E-01, 1.5700E-01, 1.5700E-01, 1.5700E-01, 1.5700E-01, 1.5700E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM + O2,O,ped,1.0000E+00,, 8.7864E+01, 8.7864E+01, 2.0318E+01, 1.7821E+01, 5.6969E+01, 5.6969E+01, 4.1130E+00, 2.0410E+00, 2.0410E+00, 1.2710E+00, 9.9600E-01, 9.9600E-01, 9.9600E-01, 9.9600E-01, 7.6200E-01, 7.6200E-01, 6.5300E-01, 6.5300E-01, 6.5300E-01, 6.5300E-01, 6.5300E-01, 6.5300E-01, 6.5300E-01, 1.1000E-02, 1.1000E-02, 1.1000E-02, 1.1000E-02, 1.1000E-02, 1.1000E-02, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, 0.0000E+00, from GITM diff --git a/src/calc_euv.cpp b/src/calc_euv.cpp index d1bb5af5..e03f045c 100644 --- a/src/calc_euv.cpp +++ b/src/calc_euv.cpp @@ -49,6 +49,9 @@ bool calc_euv(Planets planet, didWork = euv.neuvac(time, indices); else if (euvModel == "hfg") didWork = euv.solomon_hfg(time, indices); + else if (euvModel == "fism"){ + didWork = euv.get_fism(time); + } if (didWork) euv.scale_from_1au(planet, time); diff --git a/src/euv.cpp b/src/euv.cpp index e97a47a0..67865dc5 100644 --- a/src/euv.cpp +++ b/src/euv.cpp @@ -44,6 +44,11 @@ Euv::Euv() { } } + // Read in FISM data - does not need to be "slotted" + if (input.get_euv_model() == "fism") + fismData = read_fism(input.get_euv_fismfile()); + // Read in NEUVAC data - also does not need to be "slotted" + // Slot the EUVAC model coefficients: if (input.get_euv_model() == "euvac") { IsOk = slot_euv("F74113", "", euvac_f74113); @@ -167,6 +172,62 @@ bool Euv::read_file() { return DidWork; } +// ------------------------------------------------------------------------------- +// Read in FISM data. FISM files are created with srcPython/fism.py, +// and the data are read in to an index_file_output_struct. +// Inside the struct, we have time & each of the "variables" correspond to a +// FISM bin. This number of bins should match the number of bins in the EUV file +// ------------------------------------------------------------------------------- + +index_file_output_struct Euv::read_fism(std::string fism_filename) { + + std::ifstream fismfstream; + fismfstream.open(fism_filename); + std::vector> fism_file; + fism_file = read_csv(fismfstream); + + index_file_output_struct fism_contents; + + // one row per time + fism_contents.nTimes = fism_file.size(); + // first six cols are the YYYY,MM,DD,HH,mm,ss (no ms) + // the rest are the binned fism data + fism_contents.nVars = fism_file[0].size() - 6; + + // check that the user provided the correct EUV file + // The number of bins in euv file should match the number of fism bins ("nVars") + if (fism_contents.nVars != nWavelengths) { + report.error("Number of FISM wavelengths does not match the EUV file provided!"); + report.error("Either change EUV file or check your FISM file is correct."); + IsOk = false; + } + + std::vector itime(7, 0); + std::vector> values; // holds all values + std::vector values_tmp(fism_contents.nVars); // holds values in each row + + for (int iLine = 0; iLine < fism_file.size(); iLine ++) { + + itime[0] = stoi(fism_file[iLine][0]); + itime[1] = stoi(fism_file[iLine][1]); + itime[2] = stoi(fism_file[iLine][2]); + itime[3] = stoi(fism_file[iLine][3]); + itime[4] = stoi(fism_file[iLine][4]); + itime[5] = stoi(fism_file[iLine][5]); + itime[6] = 0; // 0 ms + fism_contents.times.push_back(time_int_to_real(itime)); + + for (int iVar = 0; iVar < fism_contents.nVars; iVar++) + values_tmp[iVar] = stof(fism_file[iLine][iVar + 6]); + + values.push_back(values_tmp); + } + + fism_contents.values = values; + + return fism_contents; +} + // --------------------------------------------------------------------------- // Match rows in EUV file to different types of things, such as cross // sections and spectra @@ -371,6 +432,52 @@ bool Euv::euvac(Times time, return didWork; } +// -------------------------------------------------------------------------- +// From the FISM file, interpolate the nearest 2 data to the current time +// -------------------------------------------------------------------------- + +bool Euv::get_fism(Times time) { + // This is functionally similar to get_indices, however we do not store FISM in + // the Indices class since it has variable number of bins. + + std::string function = "Euv::get_fism"; + static int iFunction = -1; + report.enter(function, iFunction); + + double time_now = time.get_current(); + bool didWork = true; + + if (fism_prev_index == 0) { + // This is probably the first time we're "running" fism. + // Make sure the file covers the entire time range of the run. + double end_time = time.get_end(); + + if (time_now < fismData.times[0] && end_time > fismData.times[-1]) { + report.error("FISM data does not cover the entire time range!"); + report.error("Please check that your FISM file is correct."); + didWork = false; + } + } + + // Get the index prior to the current time + while (fismData.times[fism_prev_index + 1] <= time_now) + fism_prev_index ++; + + // Determine time-interpolation weighting factor + precision_t dt_fism; + dt_fism = fismData.times[fism_prev_index + 1] - fismData.times[fism_prev_index]; + precision_t x = (time_now - fismData.times[fism_prev_index]) / dt_fism; + + // store the wavelength: + for (int iWave = 0; iWave < nWavelengths; iWave ++) + wavelengths_intensity_1au[iWave] = + (1.0 - x) * fismData.values[fism_prev_index][iWave] + + x * fismData.values[fism_prev_index + 1][iWave]; + + report.exit(function); + return didWork; +} + // -------------------------------------------------------------------------- // Calculate EUVAC // -------------------------------------------------------------------------- diff --git a/src/inputs.cpp b/src/inputs.cpp index 21dc05e0..74c7ed2c 100644 --- a/src/inputs.cpp +++ b/src/inputs.cpp @@ -936,13 +936,22 @@ std::string Inputs::get_diffuse_auroral_model() { } // ----------------------------------------------------------------------- -// Return the EUV model used (EUVAC only option now) +// Return the EUV model used (EUVAC, NEUVAC, FISM, etc.) // ----------------------------------------------------------------------- std::string Inputs::get_euv_model() { return mklower(check_settings_str("Euv", "Model")); } + +// ----------------------------------------------------------------------- +// Return the FISM data file +// ----------------------------------------------------------------------- + +std::string Inputs::get_euv_fismfile(){ + return get_setting_str("Euv", "fismFile"); +} + // ----------------------------------------------------------------------- // Return the heating efficiency of the neutrals for EUV // ----------------------------------------------------------------------- diff --git a/srcPython/fism.py b/srcPython/fism.py new file mode 100644 index 00000000..857845be --- /dev/null +++ b/srcPython/fism.py @@ -0,0 +1,437 @@ +#!/usr/bin/env python + +# Authors of this code: +# Daniel A. Brandt, Ph.D., Michigan Tech Research Institute, daabrand@mtu.edu +# Aaron L. Bukowski, Ph.D., University of Michigan, abukowski@umich.edu +# Aaron J. Ridley, Ph.D., University of Michigan, ridley@umich.edu + +# This file contains a suite of tools that do the following: +# 1. Download FISM2 data for a time period the user desires. +# 2. Rebin that data into the binning scheme the user desires (i.e. EUVAC-37, NEUVAC-59, or SOLOMON). +# 3. Outputs a FISM2 file with the rebinnined irradiances in the desired bins (for use by euv.cpp) + +# Top-level imports: +import argparse +import numpy as np +from datetime import datetime, timedelta +import pathlib +import os, sys +import pooch +from netCDF4 import Dataset +import scipy.integrate as integ + +# Directory management: +here = pathlib.Path(__file__).parent.resolve() +euvDir = here.parent.joinpath('share/run/UA/inputs') + +# Physical constants: +h = 6.62607015e-34 # Planck's constant in SI units of J s +c = 299792458 # Speed of light in m s^-1 + +# Helper Functions: +def getFism2(dateStart, dateEnd, source, downloadDir=None): + """ + Given a starting date and an ending date, automatically download irradiance data from LISIRD for a specific source, + including FISM2 daily or FISM2 in the Standard Bands. + :param dateStart: str + The starting date for the data in YYYY-MM-DD format. + :param dateEnd: str + The ending date for the data in YYYY-MM-DD format. + :param source: str + The type of data to be obtained. Valid inputs are: + - FISM2 (for daily averages of FISM2 data) + - FISM2S (for daily averages of FISM2 standard bands, according to Solomon and Qian 2005) + :return times: ndarray + Datetime values for each spectrum. + :return wavelengths: ndarray + Wavelength bins (bin boundaries) for the spectral data. + :return irradiance: ndarray + A 2D array where each row is a spectrum at a particular time, and the columns are wavelength bands. + """ + # Converting the input time strings to datetimes: + try: + dateStartDatetime = datetime.strptime(dateStart, "%Y-%m-%d") + dateEndDatetime = datetime.strptime(dateEnd, "%Y-%m-%d") + except: + dateStartDatetime = datetime.strptime(dateStart, "%Y%m%d") + dateEndDatetime = datetime.strptime(dateEnd, "%Y%m%d") + + # Check if the user has asked for a source that can be obtained: + validSources = ['FISM2', 'FISM2S'] + if source not in validSources: + raise ValueError("Variable 'source' must be either 'FISM2' or 'FISM2S.") + + # If the download directory is not specified, set it to the top directory that the package is in: + if downloadDir is None: + downloadDir = os.getcwd() + + # Download the most recent file for the corresponding source and read it in: + if source == 'FISM2': + url = 'https://lasp.colorado.edu/eve/data_access/eve_data/fism/daily_hr_data/daily_data.nc' + fname = 'FISM2_daily_data.nc' + urlObtain(url, loc=downloadDir, fname=fname) # hash='dbee404e1c75689b47691b8a4a733236bb66abbdc0f01b8cbd8236f69fe9d469' + datetimes, wavelengths, irradiance, uncertainties = obtainFism2(os.path.join(downloadDir, fname)) + else: + url = 'https://lasp.colorado.edu/eve/data_access/eve_data/fism/daily_bands/daily_bands.nc' + fname = 'FISM2_daily_bands.nc' + urlObtain(url, loc=downloadDir, fname=fname) # hash='27e3183f8ad6b289de191a63d3feada64c9d3f6b2973315ceda4a42c41638465' + datetimes, wavelengths, irradiance, uncertainties = obtainFism2(os.path.join(downloadDir, fname), bands=True) + + # Subset the data according to user demands: + validInds = np.where((datetimes >= dateStartDatetime) & (datetimes <= dateEndDatetime))[0] + times = datetimes[validInds] + if source == 'FISM2S': + irradiance = irradiance[-1, validInds, :] + else: + irradiance = irradiance[validInds, :] + + # Return the resulting data: + return times, wavelengths, irradiance + +def obtainFism2(myFism2File, bands=False): + """ + Load in spectrum data from a FISM2 file. + :param myFism2File: str + The location of the NETCDF4 file. + :param bands: bool + If True, loads in the data segmented into the Solomon and Qian 2005 standard bands. + :return datetimes: ndarray + An array of datetimes for each TIMED/SEE spectra. + :return wavelengths: ndarray + A one-dimensional array of wavelengths at which there are irradiance values. + :return irradiances: ndarray + A two-dimensional array of irradiance values at each time. + :return uncertainties: ndarray + A two-dimensional array of irradiance uncertainty values at each time. + """ + fism2Data = Dataset(myFism2File) + wavelengths = np.asarray(fism2Data.variables["wavelength"]) + if bands == True: # STANDARD BANDS + flux = np.asarray(fism2Data.variables["ssi"]) # photons/cm2/second + # bandwidths = np.asarray(fism2Data.variables['band_width']) + pFlux = flux * 1.0e4 # photons/m2/second + # Convert fluxes to irradiances: + irr = np.zeros_like(flux) + for i in range(flux.shape[1]): + irr[:, i] = spectralIrradiance(pFlux[:, i], wavelengths[i] * 10.0) # W/m^2 + irradiance = np.array([flux, irr]) + uncertainties = np.full_like(irradiance, fill_value=np.nan) # TODO: Replace with an estimation of uncertainty + else: # NATIVE DATA + irradiance = np.asarray(fism2Data.variables["irradiance"]) # W/m^2/nm + uncertainties = np.asarray(fism2Data.variables["uncertainty"]) + dates = fism2Data.variables["date"] + datetimes = [] + for i in range(len(dates)): + year = dates[i][:4] + day = dates[i][4:] + currentDatetime = ( + datetime(int(year), 1, 1) + + timedelta(int(day) - 1) + + timedelta(hours=12) + ) + datetimes.append(currentDatetime) + datetimes = np.asarray(datetimes) + return datetimes, wavelengths, irradiance, uncertainties + +def rebin(fism_out, saveLoc=os.getcwd(), binning_scheme='EUVAC', zero=True): + """ + Takes the output of getFism and rebins the data into whatever format the user desires. + Args: + fism_out: arraylike + The output of getFism2. Contains 4 elements: (1) datetime values for the FISM2 spectra, (2) the wavelengths + of the spectrum, (3) the actual FISM2 irradiance spectra. + saveLoc: path + Path to save data files to. Defaults to the current working directory. + binning_scheme: str + Determines the binning scheme to be used. Valid arguments include the following: + 'EUVAC' or 'Euvac' or 'euvac': Uses the 37 wavelength band scheme described in Richards, et al. 1994; doi.org/10.1029/94JA00518 + 'NEUVAC' or 'Neuvac' or 'neuvac': Uses the 59 wavelength band scheme described in Brandt and Ridley, 2024; doi.org/10.1029/2024SW004043 + 'HFG': Uses the 23 wavelength band scheme described in Solomon and Qian, 2005; https://doi.org/10.1029/2005JA011160 + 'SOLOMON' or 'Solomon' or 'solomon': Same situation as for argument 'HFG'. + NOTE: If 'HFG' or 'SOLOMON' is chosen, the values of fism_out must correspond to getFism2 being run with the + argument source='FISM2'. If this IS NOT the case, an error will be thrown. + zero: bool + Controls whether singular (bright) wavelength lines are set to a value of zero after they are extracted. + Default is True. + Returns: + fism2_file: str + The location of the rebinned FISM data. + fism2_data: arraylike + Contains 3 elements: (a) a list of datetimes for the data and (b) the rebinned FISM2 data. + """ + # Unpack the contents of fism_out: + datetimes, wavelengths, irradiance = fism_out + + # Get the native wavelength resolution of the input data: + # nativeResolution = np.concatenate((np.diff(wavelengths), np.array([np.diff(wavelengths)[-1]])), axis=0) + nativeWavelengths = wavelengths.copy() + + if binning_scheme != 'HFG' and binning_scheme != 'SOLOMON' and binning_scheme != 'Solomon' and binning_scheme != 'solomon': + if binning_scheme == 'EUVAC' or binning_scheme == 'Euvac' or binning_scheme == 'euvac': + # Grab the euv_37.csv file: + fileStr = str(euvDir.joinpath('euv.csv')) + bin_bounds = read_euv_csv_file(fileStr) + tag = '_37' + elif binning_scheme == 'NEUVAC' or binning_scheme == 'Neuvac' or binning_scheme == 'neuvac': + # Grab the euv_59.csv file: + fileStr = str(euvDir.joinpath('euv_59.csv')) + bin_bounds = read_euv_csv_file(fileStr) + tag = '_59' + else: + raise FileNotFoundError('The .csv files for specifying bin boundaries cannot be found!!') + + # Perform the rebinning! + shorts = bin_bounds['short'] / 10. + longs = bin_bounds['long'] / 10. + newWaves = 0.5 * (shorts + longs) + + # Instantiate the new data array: + if len(irradiance.shape) < 2: + fism2_data = np.zeros((1, newWaves.shape[0])) + else: + fism2_data = np.zeros((irradiance.shape[0], newWaves.shape[0])) + + # First go through all the wavelengths that are singular + myData = irradiance + for iWave, short in enumerate(shorts): + long = longs[iWave] + if (long == short): + i = np.argmin(np.abs(wavelengths - short)) + i2 = np.argmin(np.abs(nativeWavelengths - short)) + try: + fism2_data[:, iWave] = myData[:, i] * (nativeWavelengths[i2 + 1] - nativeWavelengths[i2]) + except: + fism2_data[:, iWave] = myData[i] * (nativeWavelengths[i2 + 1] - nativeWavelengths[i2]) + if zero == True: + # Zero out bin so we don't double count it. + try: + myData[:, i] = np.zeros_like(myData[:, i]) + except: + myData[i] = 0.0 + + # Then go through the ranges + for iWave, short in enumerate(shorts): + long = longs[iWave] + if (long != short): + d1 = np.abs(wavelengths - short) + iStart = np.argmin(d1) + d2 = np.abs(wavelengths - long) + iEnd = np.argmin(d2) + wave_int = 0.0 + # For wavelengths at or below 0.2 nm, just compute the sum: + if long <= 0.2: + for i in range(iStart + 1, iEnd + 1): + fism2_data[:, iWave] += myData[:, i] * \ + (wavelengths[i + 1] - wavelengths[i]) + wave_int += (wavelengths[i + 1] - wavelengths[i]) + else: + # For issues computing the sum, integrate instead: + try: + fism2_data[:, iWave] = integ.trapezoid(myData[:, iStart:iEnd], wavelengths[iStart:iEnd], axis=1) + except: + fism2_data[:, iWave] = integ.trapezoid(myData[iStart:iEnd], wavelengths[iStart:iEnd]) + + elif binning_scheme == 'HFG' or binning_scheme == 'SOLOMON' or binning_scheme == 'Solomon' or binning_scheme == 'solomon': + # Determine whether the supplied data already conforms to the Solomon and Qian binning scheme. + tag = '_solomon' + if fism_out[2].shape[1] != 23: + raise ValueError("Incorrect dimensions for element 3 of argument 'fism_out'. Dimensions must be (n,23), " + "resulting from running function 'getFism' with argument 'stanBands'=True. ") + # Should the data confirm to the proper dimensions, there is no rebinning step that needs to be done. Simply + # continue. + fism2_data = fism_out[2] + else: + # If the input irradiance data DOES NOT conform to the Solomon and Qian binning scheme, throw an error. + raise ValueError("Invalid value for argument 'binning_scheme'. Must be 'EUVAC', 'NEUVAC', 'HFG', or 'SOLOMON'.") + + # Save the rebinned data to a relative path (outside the package directory) in the form of a .txt file: + fism2_file = pathlib.Path(os.getcwd()).joinpath('fism2_file'+tag+'.txt') + saveFism(fism2_data, datetimes, fism2_file) + + return fism2_file, fism2_data + +def saveFism(data, times, filename): + """ + Takes (rebinned) FISM2 data and saves it a .txt file at a user-defined location. + Args: + data: numpy.ndarray + Irradiance data in a nxm array where the first dimension corresponds to observations (the spectrum number) + and the second dimension corresponds to wavelengths. + times: numpy.ndarray + The time values at which each spectrum is recorded. + filename: str + The desired location where the data will be saved. + Returns: + Nothing. Simply saves a file. + """ + # A helper function for working with integers: + def numStr(num): + return ',' + str(int(num)) + + # Define a helper function for opening a file to write the data, in such a way as to include parent directories if + # needed: + def safe_open_w(path): + ''' Open "path" for writing, creating any parent directories as needed. + (https://stackoverflow.com/questions/23793987/write-a-file-to-a-directory-that-doesnt-exist) + ''' + os.makedirs(os.path.dirname(path), exist_ok=True) + return open(path, 'w') + + # Open the new file and begin writing, line by line: + with safe_open_w(str(filename)) as output: + # Write the header information: + # output.write("#START\n") + # Write the irradiances themselves: + firstLine = ['%.6g' % (element) for element in data[0, :]] + firstLine_joined = ','.join(firstLine) + # The first line should always be a duplicate of the first line of data, but starting at UTC=00:00 of the first date: + output.write(str(times[0].year) + numStr(times[0].month) + numStr( + times[0].day) + ',0,0,0,' + firstLine_joined + '\n') + # The rest of the lines can be straight from the data: + for i in range(data.shape[0]): + currentLine_joined = ','.join(['%.6g' % (element) for element in data[i, :]]) + output.writelines(str(times[i].year) + numStr(times[i].month) + numStr( + times[i].day) + numStr(times[i].hour) + ',0,0,' + currentLine_joined + '\n') + # The last line should occur 12 hours from the last datapoint, but have duplicate values there: + lastLine_joined = ','.join(['%.6g' % (element) for element in data[-1, :]]) + lastTime = times[-1] + timedelta(hours=12) + output.write(str(lastTime.year) + numStr(lastTime.month) + numStr( + lastTime.day) + ',0,0,0,' + lastLine_joined + '\n') + + print('Irradiance data saved to: ') + os.system('readlink -f '+str(filename)) + return + +def spectralIrradiance(photonFlux, wavelength): + """ + Convert the photon flux to the corresponding spectral irradiance, given a specific wavelength. + Args: + photonFlux: numpy.ndarray, float, or int + Photon flux in units of photons s^-1 m^-2. For a singular wavelength, units are in photons m^-2. + wavelength: wavelength: float + A specific wavelength in Angstroms. + Returns: + irradiance: numpy.ndarray or float + The corresponding spectral irradiance in units of W/m^2/nm. + """ + photonEnergy = (h*c) / (wavelength*1e-10) # Convert the wavelength in the denominator to meters. + irradiance= photonFlux * photonEnergy + return irradiance + +def read_euv_csv_file(file): + """ + Originally written by Aaron J. Ridley, within the file 'fism2_process.py': + https://github.com/aaronjridley/EUV/blob/main/fism2_process.py + + This file reads in binning data from a CSV file that specifies bin boundaries and cross sections for either the + EUVAC model or the NEUVAC model. + Args: + file: str + The location of the .csv file to be read. + Returns: + wavelengths: numpy.ndarray + The wavelength bin boundaries for either the EUVAC model or the NEUVAC model. + """ + fpin = open(file, 'r') + + iFound = 0 + afac = [] + f74113 = [] + for line in fpin: + aline = line.split(',') + s = aline[-1].strip().split('.')[0] + if (aline[0].strip() == "Short"): + if (s.isnumeric()): + short = np.asarray(aline[5:], dtype=float) + else: + short = np.asarray(aline[5:-1], dtype=float) + iFound += 1 + if (aline[0].strip() == "Long"): + if (s.isnumeric()): + long = np.asarray(aline[5:], dtype=float) + else: + long = np.asarray(aline[5:-1], dtype=float) + if (aline[0].strip() == "F74113"): + if (s.isnumeric()): + f74113 = np.asarray(aline[5:], dtype=float) + else: + f74113 = np.asarray(aline[5:-1], dtype=float) + iFound += 1 + if (aline[0].strip() == "AFAC"): + if (s.isnumeric()): + afac = np.asarray(aline[5:], dtype=float) + else: + afac = np.asarray(aline[5:-1], dtype=float) + iFound += 1 + # Save and convert from Angstroms to nm (FISM is in nm) + wavelengths = {'short': short / 10.0, + 'long': long / 10.0, + 'afac': afac, + 'f74113': f74113} + return wavelengths + +def urlObtain(URL, loc=None, fname=None, hash=None): + """ + Helper function that uses Pooch to download files to a location specified by the user. + :param URL: str + The location of a file to be downloaded. + :param loc: str + The place the file will be downloaded. + :param fname: str + The name the file will have once it is downloaded. + :param hash: str + A known hash (checksum) of the file. Will be used to verify the download or check if an existing file needs to + be updated. + :return: + """ + if loc is None: + loc = os.getcwd() + if os.path.isfile(str(loc) + '/' + fname) is False: + fname_loc = pooch.retrieve(url=URL, known_hash=hash, fname=fname, path=loc) + else: + fname_loc = str(loc) + '/' + fname + + return fname_loc + +def get_args(): + + parser = argparse.ArgumentParser(description = 'Create FISM input data') + parser.add_argument('start', + help='Start date (format YYYYMMDD)', + type=str) + parser.add_argument('end', + help='End date (format YYYYMMDD)', + type=str) + parser.add_argument('-b', '--binning', + help="Binning scheme to use. Can be [solomon,neuvac,euvac] " + "(case insensitive)", + type=str, default="neuvac") + + args = parser.parse_args() + + return args + +# Execution (testing): +if __name__ == '__main__': + # Download some FISM2 data for the time period stated by the user. + + args = get_args() + dateStart = args.start + dateEnd = args.end + binning_scheme = args.binning + + if binning_scheme == 'HFG' or binning_scheme == 'SOLOMON' or binning_scheme == 'Solomon' or binning_scheme == 'solomon': + # SOLOMON (STAN BANDS; b23) + fism2_out_23 = getFism2(dateStart, dateEnd, 'FISM2S', downloadDir=here) + fism2_file_23, fism2_data_23 = rebin(fism2_out_23, binning_scheme=binning_scheme) + else: + fism2_out_raw = getFism2(dateStart, dateEnd, 'FISM2', downloadDir=here) + if binning_scheme == 'NEUVAC' or binning_scheme == 'Neuvac' or binning_scheme == 'neuvac': + # NEUVAC BINS (b59) + fism2_file_59, fism2_data_59 = rebin(fism2_out_raw, binning_scheme=binning_scheme, zero=True) + else: + # EUVAC BINS (b37) + fism2_file_37, fism2_data_37 = rebin(fism2_out_raw, binning_scheme=binning_scheme, zero=True) + + # Exit with a zero error code: + sys.exit(0) diff --git a/srcPython/neuvac.py b/srcPython/neuvac.py new file mode 100644 index 00000000..9134e899 --- /dev/null +++ b/srcPython/neuvac.py @@ -0,0 +1,478 @@ +#!/usr/bin/env python + +# Authors of this code: +# Daniel A. Brandt, Ph.D., Michigan Tech Research Institute, daabrand@mtu.edu + +# This file contains a suite of tools that do the following: +# 1 - Obtain F10.7 data between any dates of the user's choosing. +# 2 - Generate NEUVAC irradiances between any two dates of the user's choosing. +# 3 - Output the NEUVAC irradiances to a .csv file to be used by Aether, either in the b37 or b59 bins. + +# Top-level imports: +import argparse +import numpy as np +from datetime import datetime +import pathlib +from pathlib import Path +import pandas as pd +from scipy.interpolate import CubicSpline +import urllib.request, pickle +from fism import saveFism +import os, sys +import pooch + +# Directory management: +here = pathlib.Path(__file__).parent.resolve() +euvDir = here.parent.joinpath('share/run/UA/inputs') + +# Physical constants: +h = 6.62607015e-34 # Planck's constant in SI units of J s +c = 299792458 # Speed of light in m s^-1 + +# Global variable(s): +# Waves Table (coefficients for the old NEUVAC Model): +# Format: 0-Min, 1-Max, 2-S_1i, 3-S_Ai, 4-S_Di, 5-I_i, 6-Pi, 7-Ai +waveTable = np.array([ + [1700.00, 1750.00, 1.31491e-06, 6.71054e-06, 5.78034e-07, 0.00355128, 1.05517, 0.901612], + [1650.00, 1700.00, 5.19285e-07, 2.62376e-06, 3.08447e-07, 0.00218156, 1.06245, 0.964892], + [1600.00, 1650.00, 3.85348e-07, 1.73851e-06, 3.34911e-07, 0.00115310, 1.07246, 0.959562], + [1550.00, 1600.00, 2.96220e-07, 1.29250e-06, 2.61812e-07, 0.000814814, 1.04567, 0.967804], + [1500.00, 1550.00, 2.35326e-07, 1.21123e-06, 2.27793e-07, 0.000566574, 1.13520, 0.970257], + [1450.00, 1500.00, 1.86793e-07, 5.96399e-07, 1.48283e-07, 0.000331058, 1.01564, 0.940506], + [1400.00, 1450.00, 1.96396e-07, 5.84154e-07, 1.82438e-07, 0.000207013, 1.67546, 0.945697], + [1350.00, 1400.00, 1.04362e-07, 5.02422e-07, 1.45100e-07, 0.000153277, 1.04246, 0.992749], + [1300.00, 1350.00, 1.74403e-07, 6.32214e-07, 4.03009e-07, 0.000311075, 1.00964, 1.09381], + [1250.00, 1300.00, 7.12738e-08, 2.44220e-07, 9.56532e-08, 9.68823e-05, 1.15737, 1.01121], + [1200.00, 1250.00, 8.74335e-06, 5.02272e-05, 1.32536e-05, 0.00263307, 1.46273, 0.987493], + [1215.67, 1215.67, 6.43713e-06, 5.16823e-05, 1.11399e-05, 0.00247063, 1.26340, 0.998295], + [1150.00, 1200.00, 1.15468e-07, 2.74916e-07, 1.65125e-07, 0.000105178, 1.66887, 1.00997], + [1100.00, 1150.00, 7.71861e-08, 2.15061e-07, 1.44227e-07, 5.16157e-05, 0.971988, 1.05634], + [1050.00, 1100.00, 5.84127e-08, 3.08808e-07, 1.25160e-07, 4.65227e-05, 1.58808, 1.05327], + [1000.00, 1050.00, 2.23073e-07, 6.92710e-07, 5.19444e-07, 5.44992e-05, 0.449052, 1.10271], + [1031.91, 1031.91, 6.18723e-08, 1.21679e-07, 2.28527e-07, 3.14905e-05, 1.42684, 1.17863], + [1025.72, 1025.72, 1.61504e-07, 4.38856e-07, 2.79663e-07, 1.06365e-05, 1.09262, 1.05186], + [950.00, 1000.00, 1.70358e-07, 5.20531e-07, 3.86006e-07, 3.34989e-05, 0.491283, 1.09676], + [977.02, 977.02, 1.51857e-07, 5.60743e-07, 2.74541e-07, 6.71100e-06, 1.44918, 1.04869], + [900.00, 950.00, 7.27646e-08, 4.53511e-07, 1.91513e-07, 3.93851e-05, 1.21476, 1.06473], + [850.00, 900.00, 1.45264e-07, 2.82927e-07, 4.22856e-07, 4.83494e-05, 1.15579, 1.14948], + [800.00, 850.00, 6.69560e-08, 1.26613e-07, 1.76066e-07, 3.69687e-05, 1.14722, 1.12832], + [750.00, 800.00, 3.22816e-08, 7.81757e-08, 6.32959e-08, 4.42679e-05, 0.969748, 1.06692], + [789.36, 789.36, 1.19733e-08, 2.53334e-08, 1.58546e-08, 1.25539e-05, 1.48302, 1.00982], + [770.41, 770.41, 7.33597e-09, 2.10650e-08, 1.63125e-08, 8.88041e-06, 1.18634, 1.06584], + [765.15, 765.15, 4.85967e-09, 1.05567e-08, 5.42104e-09, 1.15262e-05, 1.17912, 1.03352], + [700.00, 750.00, 1.85139e-08, 3.63837e-08, 3.29576e-08, 1.72134e-05, 1.25328, 1.06364], + [703.36, 703.36, 5.34708e-09, 9.65120e-09, 4.54419e-09, 8.80278e-06, 1.51207, 0.972520], + [650.00, 700.00, 1.79851e-08, 6.39605e-08, 1.86000e-08, 1.41950e-05, 1.11181, 0.945801], + [600.00, 650.00, 1.52595e-07, 5.29641e-07, 1.41837e-07, 3.96165e-05, 1.00554, 0.949913], + [629.73, 629.73, 4.96048e-08, 2.46454e-07, 3.12902e-08, 1.59200e-05, 1.01611, 0.846628], + [609.76, 609.76, 2.80641e-08, 3.24530e-07, 1.81554e-08, 1.68460e-06, 0.973085, 0.793355], + [550.00, 600.00, 1.12234e-07, 6.29889e-07, 1.56092e-07, 2.79143e-05, 0.961457, 0.970150], + [584.33, 584.33, 7.91646e-08, 3.05430e-07, 5.14430e-08, 1.70372e-05, 0.844250, 0.881026], + [554.31, 554.31, 2.47485e-08, 2.68042e-07, 5.40951e-08, 1.16226e-06, 1.08699, 1.01483], + [500.00, 550.00, 1.12037e-07, 7.84515e-07, 6.32364e-08, 4.55230e-06, 1.13480, 0.816868], + [450.00, 500.00, 1.10016e-07, 3.96192e-07, 7.37101e-08, 2.62692e-05, 1.15344, 0.865234], + [465.22, 465.22, 9.60010e-09, 1.75358e-08, 6.91440e-11, 1.45142e-05, 1.62256, -0.203971], + [400.00, 450.00, 5.15555e-08, 2.89821e-07, 3.85807e-08, 1.64207e-05, 1.36652, 0.893190], + [350.00, 400.00, 3.91955e-07, 1.43942e-06, 3.16713e-07, -2.36108e-06, 1.05819, 0.910235], + [368.07, 368.07, 1.38855e-07, 7.21254e-07, 1.01814e-07, 8.71098e-07, 1.26707, 0.890513], + [300.00, 350.00, 1.35439e-06, 1.09238e-05, 8.24308e-07, 4.35250e-05, 1.22619, 0.816515], + [303.78, 303.78, 7.43959e-07, 5.94012e-06, 4.05188e-07, 9.23799e-05, 1.32976, 0.796970], + [303.31, 303.31, 5.25977e-07, 7.87164e-06, 3.07932e-07, 7.87468e-05, 0.945961, 0.759694], + [250.00, 300.00, 9.10710e-07, 3.91586e-06, 1.20177e-06, -9.64301e-06, 1.07360, 0.958369], + [284.15, 284.15, 8.67633e-07, 6.00671e-06, 3.97664e-07, -0.000107230, 1.20608, 0.773950], + [256.30, 256.30, 6.44996e-08, 4.12637e-07, 1.05193e-07, 6.61853e-06, 1.48670, 1.03265], + [200.00, 250.00, 4.83013e-07, 1.18898e-06, 8.94772e-07, 5.34779e-05, 1.04532, 1.07888], + [150.00, 200.00, 7.13305e-07, 2.47623e-06, 9.78936e-07, 0.000261230, 1.47374, 1.01156], + [100.00, 150.00, 4.03676e-08, 2.28270e-07, 4.43965e-08, 2.16162e-05, 1.09062, 0.970310], + [50.00, 100.00, 1.69769e-07, 6.93618e-07, 2.89457e-07, 2.03013e-05, 1.07887, 1.06022], + [32.00, 50.00, 1.23478e-07, 4.43644e-07, 1.75749e-07, -1.34567e-05, 1.27409, 1.01254], + [23.00, 32.00, 6.10174e-08, 2.34313e-07, 1.10591e-07, -1.22729e-05, 0.699812, 1.04841], + [16.00, 23.00, 2.23866e-07, 7.97533e-07, 3.03563e-07, -5.62012e-05, 0.706360, 0.987835], + [8.00, 16.00, 3.10773e-07, 1.22767e-06, 3.74797e-07, -8.41459e-05, 1.39529, 0.963859], + [4.00, 8.00, 1.17378e-08, 7.13970e-08, 1.38839e-08, -3.63146e-06, 0.811119, 0.920702], + [2.00, 4.00, 3.97985e-09, 4.12085e-08, 4.71914e-09, -1.86099e-06, 1.15214, 0.916686], + [1.00, 2.00, 3.52498e-09, 1.57342e-08, 4.03741e-09, -8.84488e-07, 0.951714, 0.943490] + ]) + +# Helper Functions: +def rollingAverage(myData, window_length=1, impute_edges=True, center=True): + """ + Using pandas, compute a rolling average of over 'data' using a window length of 'windowlength'. Sets the leading and + trailing windows to the values of the original data. + :param myData: arraylike + The data over which to compute the rolling average. + :param window_length: int + The size of the window over which to average. + :param impute_edges: bool + A boolean determining whether the edges will be interpolated. Default is True. + :param center: bool + A boolean determining whether the centered average will be used. + :return: rolled, arraylike + The rolling average data. + """ + myDataframe = pd.DataFrame(data=myData, columns=['Var']) + myDataframe['Rolling'] = myDataframe['Var'].rolling(window=window_length, center=center).mean() + firstValidIndex = myDataframe['Rolling'].first_valid_index() + lastValidIndex = myDataframe['Rolling'].last_valid_index() + if impute_edges == True: + # Sample x-axis: + sampleXaxis = np.linspace(0, window_length, window_length) + middleIndex = int(0.5*window_length) + # Use cubic interpolation to fill the gaps on the edges: + leadingEdgeStartingVal = myDataframe['Var'][:window_length].values[0] + leadingEndingVal = myDataframe['Rolling'][firstValidIndex] + leadingEdgeMiddleVal = np.mean([leadingEdgeStartingVal, leadingEndingVal]) + leadingSpline = CubicSpline([sampleXaxis[0], sampleXaxis[middleIndex], sampleXaxis[-1]], + [leadingEdgeStartingVal, leadingEdgeMiddleVal, leadingEndingVal]) + leadingImputedValues = leadingSpline(sampleXaxis) + + trailingEdgeStartingVal = myDataframe['Rolling'][lastValidIndex] + trailingEndingVal = myDataframe['Var'].values[-1] + trailingEdgeMiddleVal = np.mean([trailingEdgeStartingVal, trailingEndingVal]) + trailingSpline = CubicSpline([sampleXaxis[0], sampleXaxis[middleIndex], sampleXaxis[-1]], + [trailingEdgeStartingVal, trailingEdgeMiddleVal, trailingEndingVal]) + trailingImputedValues = trailingSpline(sampleXaxis) + # Ingest the imputed values: + myDataframe['Rolling'][:window_length] = leadingImputedValues + myDataframe['Rolling'][-window_length:] = trailingImputedValues + else: + myDataframe['Rolling'][:window_length] = myDataframe['Var'][:window_length] + myDataframe['Rolling'][-window_length:] = myDataframe['Var'][-window_length:] + rolled = myDataframe['Rolling'].values + return rolled + +def readCLS(filename): + """ + Load in flare-corrected, Sun-Earth distance adjusted flux values recorded by the Collecte Localisation Satellites + (CLS). + :param filename: str + The location of the data file. + :return times: list + The datetimes for each data value. + :return data: ndarray + The solar flux data for F30, F15, F10.7, F8, and F3.2. + """ + times = [] + precisionVals = [] + with open(filename, 'r') as myFile: + allLines = myFile.readlines() + data = np.zeros((len(allLines)-25, 5)) + i = 0 + j = 0 + for line in allLines: + if i >= 25: + elements = line.split() + data[j, :] = np.array([float(elements[5]), float(elements[9]), float(elements[13]), float(elements[17]), float(elements[21])]) + times.append( datetime(int(elements[0]), int(elements[1]), int(elements[2]), 12) ) + precisionVals.append( [float(elements[6]), float(elements[10]), float(elements[14]), float(elements[18]), float(elements[22])] ) + j += 1 + i += 1 + # Print the precision: + # print('Mean precision values...') + # print('F30: '+str(np.nanmean([element[0] for element in precisionVals]))+' sfu') # 6 + # print('F15: ' + str(np.nanmean([element[1] for element in precisionVals])) + ' sfu') # 8 + # print('F10.7: ' + str(np.nanmean([element[2] for element in precisionVals])) + ' sfu') # 13 + # print('F8: ' + str(np.nanmean([element[3] for element in precisionVals])) + ' sfu') # 12 + # print('F3.2: ' + str(np.nanmean([element[4] for element in precisionVals])) + ' sfu') # 11 + return times, data + +def getCLSF107(dateStart, dateEnd, truncate=True): + """ + Obtains flare-corrected F10.7 data from Collecte Localisation Satellites. The "adjusted" here means that it has been + adjusted from measurements from Earth to 1AU. (Aether/GITM need measurements at 1AU, so they can adjust to the + proper sun-planet distance, where planet can be Earth, Venus, Mars, etc.) A description of the data is + provided here: https://spaceweather.cls.fr/services/radioflux/. + Downloads the most recent measurements to a file. Reads the file and extracts the F10.7 values between two dates. + Note that if the ending date is less than or equal to the last date in the version of the file that has already + been downloaded, the file IS NOT re-downloaded, but simply parsed. Otherwise, the file is redownloaded. + :param dateStart: str + The starting date in YYYYMMDD format. + :param dateEnd: str + The ending date in YYYYMMDD format. + :param truncate: bool + Controls whether to truncate the data to exclude the most recent 81 days. Defaults is True. + """ + dateStart = dateStart[:4]+'-'+dateStart[4:6]+'-'+dateStart[6:] + dateEnd = dateEnd[:4]+'-'+dateEnd[4:6]+'-'+dateEnd[6:] + dateTimeStart = datetime.strptime(dateStart, '%Y-%m-%d') + dateTimeEnd = datetime.strptime(dateEnd, '%Y-%m-%d') + fname = euvDir.joinpath("radio_flux_adjusted_observation.txt") + if fname.exists(): + # Read in the file: + times, data = readCLS(fname) + # Check if the ending date exceeds the ending date in the file. If so, redownloading the file: + if times[-1] > dateTimeEnd: + out = urllib.request.urlretrieve( + 'ftp://ftpsedr.cls.fr/pub/previsol/solarflux/observation/radio_flux_adjusted_observation.txt', fname) + times, data = readCLS(fname) + else: + # Download the file: + out = urllib.request.urlretrieve('ftp://ftpsedr.cls.fr/pub/previsol/solarflux/observation/radio_flux_adjusted_observation.txt', fname) + times, data = readCLS(fname) + + # Compute the 81-day (centered) averaged F10.7 and 54-day averaged (: + F107 = data[:, 2] + F107A = rollingAverage(F107, window_length=81, impute_edges=True) + F107B = rollingAverage(F107, window_length=54, impute_edges=True, center=False) + print(f'\n\nlen(F107) = {len(F107)}, len(F107B) = {len(F107B)}, len(F107B) = {len(F107B)}') + # Extract the values in the desired time range: + goodInds = np.where((np.asarray(times) >= dateTimeStart) & (np.asarray(times) <= dateTimeEnd))[0] + # Truncation: + if truncate and len(goodInds) >= 2*81: + goodInds = goodInds[:-81] + return np.asarray(times)[goodInds], np.asarray(F107)[goodInds], np.asarray(F107A)[goodInds], np.asarray(F107B)[goodInds] + +def mycorrelate2d(df, normalized=False): + """ + Compute the correlation matrix from 2D data, where each row is cross correlated with the others. + This function handles NaN values by ignoring them. + :param df: ndarray + A 2D array of dimensions n x m. + :param normalized: bool + Determines whether the resulting correlation matrix is normalized. Default is False. + :returns ccm: ndarray + The [normalized] cross-correlation matrix. + Source: https://stackoverflow.com/questions/54292947/basics-of-normalizing-cross-correlation-with-a-view-to-comparing-signals + """ + # Initialize cross correlation matrix with zeros + ccm = np.zeros((df.shape[1], df.shape[1])) + # Fill in each entry of the matrix one-by-one: + for i in range(df.shape[1]): + outer_row = df[:, i] + for j in range(df.shape[1]): + inner_row = df[:, j] + goodInds = np.logical_and(~np.isnan(outer_row), ~np.isnan(inner_row)) + if (not normalized): + x = np.correlate(outer_row[goodInds], inner_row[goodInds]) + else: + x = get_cc(outer_row[goodInds], inner_row[goodInds]) + # a = (inner_row - np.mean(inner_row)) / (np.std(inner_row) * len(inner_row)) + # b = (outer_row - np.mean(outer_row)) / (np.std(outer_row) ) + # x = np.correlate(a, b) + ccm[i, j] = x + return ccm + +def get_cc(array1, array2, normalize=True): + """ + Compute the cross-correlation of two 1D arrays of the same length. + :param array1: ndarray + A 1D array of length n. + :param array2: ndarray + A 1D array of length n. + :return c: float + The normalized correlation of the two arrays. + """ + if normalize: + a = (array1 - np.mean(array1)) / (np.std(array1) * len(array1)) + b = (array2 - np.mean(array2)) / (np.std(array2)) + c = np.correlate(a, b) + else: + c = np.correlate(array1, array2) + return c + +def loadPickle(pickleFilename): + """ + Given the name of a (pre-existing) pickle file, load its contents. + :param: pickleFilename, str + A string with the location/name of the filename. + :return: var + The loaded data. + """ + with open(pickleFilename, 'rb') as pickleFile: + var = pickle.load(pickleFile) + return var + +# CORE NEUVAC FUNCTIONS: +def irrFunc(F107input, A, B, C, D, E, F): + F107, F107A = F107input + return A * (F107 ** B) + C * (F107A ** D) + E * (F107A - F107) + F + +def neuvacEUV(f107, f107b, bands=None, tableFile=None, statsFiles=None): + """ + Use a parametric model to compute solar flux in the 59 conventional wavelength bands used by Aether/GITM. Capable + of returning perturbed irradiance values that are perturbed according to the variations of in the intensity of each + bin + :param f107: ndarray + F10.7 values. + :param f107b: ndarray + 81-day center-averaged F10.7 values; must be the same length as f107. + :param bands: str + If None or 'NEUVAC', returns irradiances in the GITM Bands. If 'EUVAC', returns them in the 37 bands used by + EUVAC. If 'SOLOMON', returns them in the 22 bands used by Solomon and Qian. + :param tableFile: str + Corresponds to the .txt file holding the NEUVAC coefficients most recently-generated by fitNeuvac.py. If + not given, simply uses the table file corresponding to the selected bin structure. Default is None. + :param statsFiles: Bool + Determines whether data for uncertainty quantification is exploited. Involves usage of a list containing + 2 elements where the first element is a file containing the 59x59 correlation matrix and the second + element is a file containing the 1x59 standard deviation values for NEUVAC. NOT REQUIRED. + :return euvIrradiance: ndarray + A nxm ndarray where n is the number of EUV irradiance values and m is the number of wavelength bands. + :return perturbedEuvIrradiance: ndarray + A nxm ndarray where n is the number of EUV irradiance values perturbed due to inherent uncertainty and m is the + number of wavelength bands. + :return savedPerts: ndarray + A nxm ndarray of the perturbations (time series of the NEUVAC+Perturbation - NEUVAC) + :return cc2: ndarray + A mxm ndarray of the correlation matrix between each wavelength's time-series of the NEUVAC+Perturbation - + NEUVAC. + """ + if type(f107) != np.ndarray: + f107 = np.asarray([f107]) + f107b = np.asarray([f107b]) + if bands == 'SOLOMON': + solarFlux = np.zeros((1, 22)) + else: + solarFlux = np.zeros((1, waveTable.shape[0])) + else: + if bands == 'SOLOMON': + solarFlux = np.zeros((len(f107), 22)) + else: + solarFlux = np.zeros((len(f107), waveTable.shape[0])) + euvIrradiance = np.zeros_like(solarFlux) + perturbedEuvIrradiance = np.zeros_like(solarFlux) + # Gather the model parameters: + if tableFile is None: + if bands == 'SOLOMON': + tableFile = euvDir.joinpath('neuvac_table_stan_bands.txt') #'../data/neuvac_table_stan_bands.txt' + else: + tableFile = euvDir.joinpath('neuvac_table.txt') #'../data/neuvac_table.txt' + neuvacTable = [] + with open(here.parent.joinpath(tableFile)) as neuvacFile: # open(tableFile) + contents = neuvacFile.readlines() + i = 0 + for line in contents: + if i > 17: + neuvacTable.append([float(element) for element in line.split(' ')]) + i+=1 + neuvacTable = np.asarray(neuvacTable) + + # If no stats file is provided, simply return the base model output (using the required table file): + if not statsFiles: + # Loop across the F10.7 (and F10.7A) values: + for i in range(len(f107)): + k = 0 + for j in (range(solarFlux.shape[1])): + irrRes = irrFunc([f107[i], f107b[i]], *neuvacTable[j, 2:]) + if irrRes < 0: + irrRes = 0 + euvIrradiance[i, k] = irrRes + else: + euvIrradiance[i, k] = irrRes + k += 1 + if bands == 'EUVAC': # Returns values ONLY for those corresponding to the wavelengths used by EUVAC + return euvIrradiance[:, 7:44], None, None, None + else: + return euvIrradiance, None, None, None + else: + # Include statistical data for calculating uncertainties via perturbations: + if bands == 'NEUVAC': + statsFiles = [euvDir.joinpath('corMat.pkl'), euvDir.joinpath('sigma_NEUVAC.pkl')] + elif bands == 'EUVAC': + statsFiles = [euvDir.joinpath('corMatEUVAC.pkl'), euvDir.joinpath('sigma_EUVAC.pkl')] + else: + statsFiles = [euvDir.joinpath('corMatStanBands.pkl'), euvDir.joinpath('sigma_NEUVAC_StanBands.pkl')] + corMatFile = statsFiles[0] + corMat = loadPickle(corMatFile) + sigmaFile = statsFiles[1] + STDNeuvacResids = loadPickle(sigmaFile) + # Loop across the F10.7 (and F10.7A) values: + nTimes = len(f107) + nWaves = solarFlux.shape[1] + savedPerts = np.zeros((nTimes, nWaves)) + for i in range(len(f107)): + # Loop across the wavelengths (59 conventional wavelengths): + k = 0 + P_n = [] + for j in (range(solarFlux.shape[1])): + # Percentage perturbation: + P_j = np.random.normal(0, 1.0) + P_n.append(P_j) + P_1 = P_n[0] + # Normalized Correlated Perturbation: + if bands == 'SOLOMON': + if j < 5: + C_j1 = corMat[0, j] # 3 # Only consider correlation with the third wavelength bin of the SOLOMON bins! + else: + C_j1 = corMat[5, j] # Only consider correlation with the fifth wavelength bin of the SOLOMON bins! + else: + if j < 7: + # Only consider correlation with the third wavelength bin (of the NEUVAC bins!) when bands are below 8. + C_j1 = corMat[0, j] # 2 + else: + # Only consider correlation with the first wavelength bin (of the EUVAC bins!) when bands are above 8. + C_j1 = corMat[7, j] + N_j = C_j1 * P_1 + (1.0 - C_j1) * P_j + # Actual Normalized Correlated Perturbation: + A_j = STDNeuvacResids[j] * N_j + irrRes = irrFunc([f107[i], f107b[i]], *neuvacTable[j, 2:]) + if irrRes < 0: + irrRes = 0 + euvIrradiance[i, k] = irrRes + if irrRes + A_j < 0: + perturbedEuvIrradiance[i, k] = 0 + else: + perturbedEuvIrradiance[i, k] = irrRes + A_j + else: + euvIrradiance[i, k] = irrRes + perturbedEuvIrradiance[i, k] = irrRes + A_j + savedPerts[i, j] = A_j + k += 1 + + # Generate a correlation matrix of the perturbations (to compare to the input correlation matrix as a sanity check): + cc2 = mycorrelate2d(savedPerts, normalized=True) + + if bands == 'EUVAC': # Returns values ONLY for those corresponding to the wavelengths used by EUVAC + return euvIrradiance[:, 7:44], perturbedEuvIrradiance[:, 7:44], savedPerts, cc2 + else: + return euvIrradiance, perturbedEuvIrradiance, savedPerts, cc2 + +# ----------------------------------------------------------------------------------------------------------------------------------------- +# Argument Parsing Function: +def get_args(): + + parser = argparse.ArgumentParser(description = 'Create NEUVAC input data') + parser.add_argument('start', + help='Start date (format YYYYMMDD)', + type=str) + parser.add_argument('end', + help='End date (format YYYYMMDD)', + type=str) + parser.add_argument('-b', '--binning', + help="Binning scheme to use. Can be [solomon,neuvac,euvac] " + "(case insensitive)", + type=str, default="neuvac") + args = parser.parse_args() + + return args +# ----------------------------------------------------------------------------------------------------------------------------------------- + +# Example Execution: + +# python neuvac.py 20110319 20110321 -b euvac + +args = get_args() +dateStart = args.start +dateEnd = args.end +binning_scheme = args.binning + +# Load F10.7 data (from Collecte Localisation Satellites): +times, F107, F107A, F107B = getCLSF107(dateStart, dateEnd) + +# Generate NEUVAC Irradiance from that F10.7 data: +if binning_scheme == 'HFG' or binning_scheme == 'SOLOMON' or binning_scheme == 'Solomon' or binning_scheme == 'solomon': + # SOLOMON (STAN BANDS; b23) + irradiance, _, _, _ = neuvacEUV(F107, F107B, bands='SOLOMON') +else: + if binning_scheme == 'NEUVAC' or binning_scheme == 'Neuvac' or binning_scheme == 'neuvac': + # NEUVAC BINS (b59) + irradiance, _, _, _ = neuvacEUV(F107, F107B, bands='NEUVAC') + else: + # EUVAC BINS (b37) + irradiance, _, _, _ = neuvacEUV(F107, F107B, bands='EUVAC') + +# Save the NEUVAC Irradiance to a file that Aether can use: +tag = str(irradiance.shape[1]) +fname = os.getcwd() + '/neuvac_file_'+tag+'.txt' +saveFism(irradiance, times, fname) +