PyLAP changes to work with sigmond.

.gitignore

@@ -1,9 +1,13 @@
 *.pyc
-
+*.so
+build/
+venv/
+__pycache__/
 
 Ionosphere/__pycache__/*
 
 Examples/results/*
+iri2020_test_output.png
 *.nc
 .nc
 *.csv

Examples/SAMI3_ray_test1.py

@@ -9,7 +9,6 @@
 from Plotting import plot_ray_iono_slice as plot_iono
 import matplotlib.pyplot as plt
 import datetime as dt
-import ipdb
 plt.switch_backend('tkagg')
 
 # Constants

Examples/run_raytraces.py

@@ -11,7 +11,6 @@
 import matplotlib.pyplot as plt
 import datetime as dt
 import pandas as pnd
-import ipdb
 import os
 from geographiclib.geodesic import Geodesic
 geod = Geodesic.WGS84

Examples/test_iri2020.py

@@ -0,0 +1,89 @@
+#!/usr/bin/env python3
+"""
+Test IRI-2020 ionosphere model (works without PHaRLAP)
+
+This script demonstrates the IRI-2020 integration and can be run
+before PHaRLAP is installed to verify the ionosphere model works.
+"""
+
+import datetime
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Import iri2020 package
+import iri2020
+
+def main():
+    # Test parameters
+    lat = 40.0      # Latitude (degrees)
+    lon = -75.0     # Longitude (degrees)
+    dt = datetime.datetime(2024, 6, 21, 12, 0)  # Summer solstice, noon UTC
+
+    # Height range (km) - IRI2020 expects [min, max, step]
+    alt_min = 100
+    alt_max = 600
+    alt_step = 5
+    heights = [alt_min, alt_max, alt_step]
+
+    print(f"Running IRI-2020 for:")
+    print(f"  Location: {lat}°N, {lon}°E")
+    print(f"  Time: {dt}")
+    print(f"  Heights: {alt_min} - {alt_max} km (step: {alt_step} km)")
+    print()
+
+    # Call IRI-2020
+    print("Calling IRI-2020...")
+    result = iri2020.IRI(dt, heights, lat, lon)
+
+    # Reconstruct height array for plotting
+    height_arr = np.arange(alt_min, alt_max + alt_step, alt_step)
+
+    # Extract electron density
+    ne = result['ne'].values  # electrons/m^3
+
+    # Find F2 peak
+    f2_idx = np.nanargmax(ne)
+    f2_height = height_arr[f2_idx]
+    f2_density = ne[f2_idx]
+
+    # Calculate plasma frequency at F2 peak
+    # fp = sqrt(ne * e^2 / (4 * pi^2 * epsilon_0 * m_e))
+    # Simplified: fp (MHz) = 9e-6 * sqrt(ne)
+    fp_f2 = 9e-6 * np.sqrt(f2_density)
+
+    print(f"Results:")
+    print(f"  F2 peak height: {f2_height:.1f} km")
+    print(f"  F2 peak density: {f2_density:.2e} electrons/m³")
+    print(f"  F2 critical frequency (foF2): {fp_f2:.2f} MHz")
+    print()
+
+    # Plot
+    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
+
+    # Electron density profile
+    ax1.plot(ne / 1e11, height_arr, 'b-', linewidth=2)
+    ax1.axhline(y=f2_height, color='r', linestyle='--', label=f'F2 peak: {f2_height:.0f} km')
+    ax1.set_xlabel('Electron Density (×10¹¹ m⁻³)')
+    ax1.set_ylabel('Height (km)')
+    ax1.set_title(f'IRI-2020 Electron Density Profile\n{lat}°N, {lon}°E, {dt}')
+    ax1.legend()
+    ax1.grid(True, alpha=0.3)
+
+    # Temperature profiles
+    Te = result['Te'].values
+    Ti = result['Ti'].values
+    ax2.plot(Te, height_arr, 'r-', linewidth=2, label='Electron Temp (Te)')
+    ax2.plot(Ti, height_arr, 'b-', linewidth=2, label='Ion Temp (Ti)')
+    ax2.set_xlabel('Temperature (K)')
+    ax2.set_ylabel('Height (km)')
+    ax2.set_title('IRI-2020 Temperature Profiles')
+    ax2.legend()
+    ax2.grid(True, alpha=0.3)
+
+    plt.tight_layout()
+    plt.savefig('iri2020_test_output.png', dpi=150)
+    print("Plot saved to: iri2020_test_output.png")
+    plt.show()
+
+if __name__ == '__main__':
+    main()

Ionosphere/gen_SAMI3_iono_grid_2d.py

@@ -1,7 +1,6 @@
 import netCDF4 as nc
 import scipy as sp
 import numpy as np
-import ipdb
 import datetime as dt
 import pandas as pd
 import tqdm

Ionosphere/gen_iono_grid_2d.py

@@ -167,6 +167,8 @@
 from pylap.iri2016 import iri2016
 from pylap.iri2012 import iri2012
 from pylap.iri2007 import iri2007
+import iri2020 as iri2020_pkg
+import datetime
 #
 #function [iono_pf_grid,iono_pf_grid_5,collision_freq,irreg,iono_te_grid] = ...
 #    gen_iono_grid_2d(origin_lat, origin_lon, R12, UT, azim, ...
@@ -194,14 +196,11 @@ def gen_iono_grid_2d(origin_lat, origin_lon, R12, UT, azim,
     # defalut value
 
 
-    if profile_type.lower() != 'chapman_fllhc' and \
-    	profile_type.lower() != 'chapman' and \
-    	profile_type.lower() != 'iri' and \
-      profile_type.lower() != 'iri2007'     and \
-    	profile_type.lower() !=  'iri2012' and \
-    	profile_type.lower() !=  'iri2016'  and \
-    	profile_type.lower() !=  'firi':
-        print('invalid profile type')
+    valid_profile_types = ['chapman_fllhc', 'chapman', 'iri', 'iri2007', 
+                             'iri2012', 'iri2016', 'iri2020', 'firi']
+    if profile_type.lower() not in valid_profile_types:
+        print('invalid profile type: {}. Valid types: {}'.format(
+              profile_type, valid_profile_types))
         sys.exit('gen_iono_grid_2d')
 #*
     fllhc_flag = 0
@@ -380,7 +379,7 @@ def gen_iono_profile(lat, lon, num_heights, start_height, height_inc,
       #M%%%%%%%%%%%%%%%%%
       #MThis is IRI2012 %
       #M%%%%%%%%%%%%%%%%%
-  elif profile_type.lower == 'iri2012':
+  elif profile_type.lower() == 'iri2012':
 
   #        #M call IRI 2012	   
     iono, iono_extra = iri2012.iri2012(lat, lon, R12, UT, start_height, 
@@ -417,7 +416,7 @@ def gen_iono_profile(lat, lon, num_heights, start_height, height_inc,
       #M%%%%%%%%%%%%%%%%%
       #MThis is IRI2007 %
       #M%%%%%%%%%%%%%%%%%
-  elif profile_type.lower == 'iri2007':
+  elif profile_type.lower() == 'iri2007':
     #M IRI2007 only returns 100 values for electron density with height - so
     #M determine the number of  multiple calls required.
     max_iri_numhts = 100
@@ -471,6 +470,48 @@ def gen_iono_profile(lat, lon, num_heights, start_height, height_inc,
       iono_ti_prof[idx] = ion_temp
       iono_ti_prof[iono_ti_prof == -1] = np.nan
       iono_te_prof[iono_te_prof == -1] = np.nan
+
+      #M%%%%%%%%%%%%%%%%%
+      #MThis is IRI2020 %
+      #M%%%%%%%%%%%%%%%%%
+  elif profile_type.lower() == 'iri2020':
+    # Use the iri2020 Python package (space-physics/iri2020)
+    # Convert UT array [year, month, day, hour, minute] to datetime
+    dt_time = datetime.datetime(UT[0], UT[1], UT[2], UT[3], UT[4])
+
+    # Generate height array
+    height_arr = np.arange(num_heights) * height_inc + start_height
+
+    # Call IRI2020 - returns xarray Dataset
+    iri_result = iri2020_pkg.IRI(dt_time, height_arr, lat, lon)
+
+    # Extract electron density (m^-3)
+    elec_dens = iri_result['ne'].values
+    elec_dens[np.isnan(elec_dens)] = 0
+    elec_dens[elec_dens < 0] = 0
+
+    # Extract temperatures
+    iono_te_prof = iri_result['Te'].values
+    iono_ti_prof = iri_result['Ti'].values
+    iono_te_prof[iono_te_prof < 0] = np.nan
+    iono_ti_prof[iono_ti_prof < 0] = np.nan
+
+    # Calculate plasma frequency (MHz)
+    iono_pf_prof = np.sqrt(pfsq_conv * elec_dens)
+
+    if doppler_flag:
+      # UT 5 minutes later
+      dt_time_5 = dt_time + datetime.timedelta(minutes=5)
+      iri_result_5 = iri2020_pkg.IRI(dt_time_5, height_arr, lat, lon)
+      elec_dens5 = iri_result_5['ne'].values
+      elec_dens5[np.isnan(elec_dens5)] = 0
+      elec_dens5[elec_dens5 < 0] = 0
+      iono_pf_prof5 = np.sqrt(pfsq_conv * elec_dens5)
+    else:
+      iono_pf_prof5 = iono_pf_prof.copy()
+
+    # iono_extra placeholder - IRI2020 package returns different structure
+    iono_extra = None
 
   print('leave gen_iono_profile')
 

Ionosphere/gen_iono_grid_3d.py

@@ -154,6 +154,8 @@
 from pylap.iri2016 import iri2016
 from pylap.iri2012 import iri2012
 from pylap.iri2007 import iri2007
+import iri2020 as iri2020_pkg
+import datetime
 
 def gen_iono_grid_3d(UT, R12, iono_grid_parms,
                      geomag_grid_parms, doppler_flag, 
@@ -167,15 +169,12 @@ def gen_iono_grid_3d(UT, R12, iono_grid_parms,
     else:
         iri_options = {}
 
-    if profile_type.lower() != 'chapman_fllhc' and \
-    	profile_type.lower() != 'chapman' and \
-    	profile_type.lower() != 'iri' and \
-      profile_type.lower() != 'iri2007'     and \
-    	profile_type.lower() !=  'iri2012' and \
-    	profile_type.lower() !=  'iri2016'  and \
-    	profile_type.lower() !=  'firi':
-        print('invalid profile type')
-        sys.exit('gen_iono_grid_2d')
+    valid_profile_types = ['chapman_fllhc', 'chapman', 'iri', 'iri2007', 
+                             'iri2012', 'iri2016', 'iri2020', 'firi']
+    if profile_type.lower() not in valid_profile_types:
+        print('invalid profile type: {}. Valid types: {}'.format(
+              profile_type, valid_profile_types))
+        sys.exit('gen_iono_grid_3d')
 #*
     fllhc_flag = 0
     if profile_type.lower() == 'chapman_fllhc':
@@ -365,7 +364,7 @@ def gen_iono_subgrid(lat, lon_min, lon_inc, lon_max, ht_min, ht_inc,
 #     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 #     % This is IRI2016 with FIRI option on %
 #     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-        elif profile_type.lower == 'firi':
+        elif profile_type.lower() == 'firi':
 
             num_heights = round((ht_max - ht_min)/ht_inc) + 1
 
@@ -520,5 +519,57 @@ def gen_iono_subgrid(lat, lon_min, lon_inc, lon_max, ht_min, ht_inc,
                 collision_freq_subgrid[lon_idx,:] = \
                 eff_coll_freq.eff_coll_freq(T_e, T_ion, elec_dens, neutral_dens)
 
+#     %%%%%%%%%%%%%%%%%%%
+#     % This is IRI2020 %
+#     %%%%%%%%%%%%%%%%%%%
+        elif profile_type.lower() == 'iri2020':
+            # Use the iri2020 Python package (space-physics/iri2020)
+            num_heights = round((ht_max - ht_min)/ht_inc) + 1
+            height_arr = np.arange(num_heights) * ht_inc + ht_min
+
+            # Convert UT array [year, month, day, hour, minute] to datetime
+            dt_time = datetime.datetime(UT[0], UT[1], UT[2], UT[3], UT[4])
+
+            # Call IRI2020 - returns xarray Dataset
+            iri_result = iri2020_pkg.IRI(dt_time, height_arr, lat, lon)
+
+            # Extract electron density (m^-3)
+            elec_dens = iri_result['ne'].values
+            elec_dens[np.isnan(elec_dens)] = 0
+            elec_dens[elec_dens < 0] = 0
+
+            # Calculate plasma frequency (MHz)
+            iono_pf_subgrid[lon_idx] = np.sqrt(elec_dens * pfsq_conv)
+
+            if doppler_flag:
+                dt_time_5 = dt_time + datetime.timedelta(minutes=5)
+                iri_result_5 = iri2020_pkg.IRI(dt_time_5, height_arr, lat, lon)
+                elec_dens5 = iri_result_5['ne'].values
+                elec_dens5[np.isnan(elec_dens5)] = 0
+                elec_dens5[elec_dens5 < 0] = 0
+                iono_pf_subgrid_5[lon_idx] = np.sqrt(elec_dens5 * pfsq_conv)
+            else:
+                iono_pf_subgrid_5[lon_idx] = np.sqrt(elec_dens * pfsq_conv)
+
+            # Extract temperatures
+            T_e = iri_result['Te'].values
+            T_e[T_e < 0] = np.nan
+
+            T_ion = iri_result['Ti'].values
+            T_ion[T_ion < 0] = np.nan
+
+            # Neutral densities
+            lat_arr = lat * np.ones_like(height_arr)
+            lon_arr = lon * np.ones_like(height_arr)
+            if R12 == -1:
+                [neutral_dens, temp] = nrlmsise00(lat_arr, lon_arr, height_arr, UT)
+            else:
+                f107 = 63.75 + R12*(0.728 + R12*0.00089)
+                [neutral_dens, temp] = nrlmsise00(lat_arr, lon_arr, height_arr, UT, f107, f107, 4)
+
+            # Calculate collision frequency
+            collision_freq_subgrid[lon_idx,:] = \
+            eff_coll_freq.eff_coll_freq(T_e, T_ion, elec_dens, neutral_dens)
+
 
     return iono_pf_subgrid, iono_pf_subgrid_5, collision_freq_subgrid

Plotting/Plot_map.py

@@ -7,7 +7,6 @@
 import os
 import netCDF4 as nc
 import datetime as dt
-import ipdb
 
 
 def calculate_scale(data,stddevs=2.,lim='auto'):

Plotting/plot_ray_iono_slice.py

@@ -74,7 +74,6 @@
 import sys
 import platform
 from qtpy.QtWidgets import QApplication
-import ipdb
 
 
 def plot_ray_iono_slice(iono_grid, start_range, end_range, range_inc,