diff --git a/.gitignore b/.gitignore
index a52a30c..27a8516 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,9 +1,13 @@
 *.pyc
-
+*.so
+build/
+venv/
+__pycache__/
 
 Ionosphere/__pycache__/*
 
 Examples/results/*
+iri2020_test_output.png
 *.nc
 .nc
 *.csv
diff --git a/Examples/SAMI3_ray_test1.py b/Examples/SAMI3_ray_test1.py
index 15e89d2..3b2fca1 100644
--- a/Examples/SAMI3_ray_test1.py
+++ b/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
diff --git a/Examples/run_raytraces.py b/Examples/run_raytraces.py
index 55a3f89..709b60b 100644
--- a/Examples/run_raytraces.py
+++ b/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
diff --git a/Examples/test_iri2020.py b/Examples/test_iri2020.py
new file mode 100644
index 0000000..bfa6324
--- /dev/null
+++ b/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()
diff --git a/Ionosphere/gen_SAMI3_iono_grid_2d.py b/Ionosphere/gen_SAMI3_iono_grid_2d.py
index 18531ab..ab1fa7c 100644
--- a/Ionosphere/gen_SAMI3_iono_grid_2d.py
+++ b/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
diff --git a/Ionosphere/gen_iono_grid_2d.py b/Ionosphere/gen_iono_grid_2d.py
index 1924520..304e9ef 100644
--- a/Ionosphere/gen_iono_grid_2d.py
+++ b/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')
 
diff --git a/Ionosphere/gen_iono_grid_3d.py b/Ionosphere/gen_iono_grid_3d.py
index 95d73ad..fc30604 100644
--- a/Ionosphere/gen_iono_grid_3d.py
+++ b/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
diff --git a/Plotting/Plot_map.py b/Plotting/Plot_map.py
index 36240ae..0da9ace 100644
--- a/Plotting/Plot_map.py
+++ b/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'):
diff --git a/Plotting/plot_ray_iono_slice.py b/Plotting/plot_ray_iono_slice.py
index be5fb6d..546f089 100644
--- a/Plotting/plot_ray_iono_slice.py
+++ b/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,
diff --git a/README.md b/README.md
index d708213..a6e2442 100644
--- a/README.md
+++ b/README.md
@@ -1,76 +1,166 @@
-# Experimental Software Warning
+# PyLap
 
-PyLap is based on the scientific software PHaRLAP and is not meant for operational use as is stated on PHaRLAP's website. any and all risk falls to you if you are using this software.
+A Python 3 wrapper for the PHaRLAP ionospheric raytracer.
 
-## Linux Install Instructions ##
+> **Warning:** PyLap is based on the scientific software PHaRLAP and is not meant for operational use as stated on PHaRLAP's website. Any and all risk falls to you if you are using this software.
 
-*Note* This install is only available for Ubuntu linux systems with X86 CPU's *Note* 
+## About this fork
 
-*Note*  It is also possible to install PyLap using WSL2 with Ubuntu for windows machines running windows 10 and windows 11 *Note* 
+This is a patched fork of [HamSCI/PyLap](https://github.com/HamSCI/PyLap) maintained for use with [hf-timestd](https://github.com/mijahauan/hf-timestd) and other HF-propagation projects. The upstream version will not build or run correctly against current PHaRLAP (4.7.4) without these patches:
 
-Download Steps
+- **PHaRLAP 4.7.4 support** — upstream targets 4.5.x; this fork's `setup.py` links against 4.7.4 static libraries and gracefully skips legacy IRI modules (`iri2007`, `iri2012`) when their `.a` files are absent.
+- **Cross-platform builds** — adds macOS arm64 and macOS x86_64 alongside Linux x86_64 (upstream is Linux-only).
+- **Unified gfortran build** — Intel Fortran redistributable is no longer required; Linux and macOS both build against gfortran, eliminating a significant install-time hurdle.
+- **GCC 14+ compatibility** — adds `-Wno-error=incompatible-pointer-types` so builds succeed on Debian 13/trixie and other distros that ship GCC 14 (which promotes this warning to an error).
+- **Multi-hop stride fix in `raytrace_2d.c`** — the `ray_data` C-array hop stride was `num_rays × 9`; corrected to `num_rays × 24` so fields past hop 0 are read from the right offsets under PHaRLAP 4.7.4.
+- **`iri2016` module now calls IRI-2020** — PHaRLAP 4.7.4 replaced IRI-2016 internals with IRI-2020; the fork's wrapper was updated to match.
+- **Fortran SAVE-variable segfault workaround (caller-side note)** — repeated `raytrace_2d` calls can crash due to persistent Fortran state; make a single call with `nhops=max_hops` and iterate hops in Python.
+- **Ne unit convention (caller-side note)** — IRI-2020 returns electron density in m⁻³ but `raytrace_2d` expects cm⁻³; scale by 10⁻⁶ before passing the grid.
 
-1. Create a folder on the home directory where you will place all of the downloaded resources.
+## Requirements
 
-2. Download PHaRLAP toolbox for matlab and unzip to the directory that you just created https://www.dst.defence.gov.au/our-technologies/pharlap-provision-high-frequency-raytracing-laboratory-propagation-studies .
+- **OS:** Ubuntu/Debian Linux (x86_64), WSL2 on Windows 10/11, macOS (arm64 or x86_64)
+- **Python:** 3.8+
+- **PHaRLAP:** 4.7.4 required for raytracing (request access from DST Australia — see below)
+- **gfortran:** Required to build PyLap and the IRI-2020 Python package — `apt install gfortran` on Linux, `brew install gcc` on macOS. Intel Fortran is **not** required (removed in this fork).
 
-3. Download Pylap from github in the same directory that you just created. 
+## Quick Start
 
-    *Note* This works for either cloning the github repository or downloading it as a zip file *Note*
+### 1. Obtain PHaRLAP
 
-4. Download Redistributable Libraries for Intel® C++ and Fortran 2020. Which is required because the original fortran code was compiled using the Intel fortran compiler. The one used originally is available at the following download link: https://www.intel.com/content/www/us/en/developer/articles/tool/redistributable-libraries-for-intel-c-and-fortran-2020-compilers-for-linux.html
+PHaRLAP is **not** redistributed with this repo — you must obtain it separately.
 
-5. cd within a terminal window into the intel libraries folder and run install.sh, follow the prompt until install complete. when promted with where to install this, the default should be the home drive (ex. [/home/UserName]). The default location is fine to use for the installation.
+**PHaRLAP 4.7.4** (required for raytracing):
+- Request access at https://www.dst.defence.gov.au/our-technologies/pharlap-provision-high-frequency-raytracing-laboratory-propagation-studies
+- DST will email you a download link after reviewing your request (typically 1–3 business days)
+- License restricts use to research; see PHaRLAP's own license terms
+- After extracting, verify `lib/` contains `linux/`, `maca/`, `maci/` subdirectories (these hold the static libraries this fork links against)
 
+Suggested layout:
 
-File Directory model
+```
+~/pylap_project/
+├── PyLap/            # This repository
+└── pharlap_4.7.4/    # PHaRLAP toolbox (from DST)
+```
 
-├── PyLap Project folder  
-   &emsp;├── PyLap\
-   &emsp;├── PHARLAP\
-   &emsp;└── l_comp_lib_2020.4.304_comp.for_redist\
+### 2. Run Setup Script
 
+```bash
+cd ~/pylap_project/PyLap
+. ./setup.sh
+# Enter: ~/pylap_project
+```
 
-## Setup with script file
+The setup script will:
+- Create a Python virtual environment at `PyLap/venv/`
+- Install system dependencies (requires sudo)
+- Install Python packages from `requirements.txt`
+- Build and install PyLap
+- Add a `pylap-activate` alias to your `.bashrc`
 
-1. cd into the pylap project directory.
+### 3. Activate and Run
 
-2. Run the setup.sh script using the command “. ./setup.sh”. running this script will promt you to enter the filepath of the folder that all of your project is installed.
+```bash
+# Activate the virtual environment
+source venv/bin/activate
+# Or use the alias:
+pylap-activate
 
-3. everything should be setup correctly! run the example files that are in the examples folder to make sure everything is setup correctly! if for some reason the example files do not work check the bashrc file in your home directory by using the command "nano .bashrc"(must be in the home directory).
-  
-    *Note* This is a one time setup and does not need to be run again unless a new install is made *Note*
+# Test IRI-2020 (works without PHaRLAP)
+python3 Examples/test_iri2020.py
 
+# Test raytracing (requires PHaRLAP)
+python3 Examples/ray_test1.py
+```
 
+> **Note:** The raytracing examples (`ray_test1.py`, etc.) require PHaRLAP to be installed. Use `test_iri2020.py` to verify your installation before PHaRLAP is available.
 
-## Manual setup
+### 4. Verify the Build
 
-5. export PHARLAP_HOME="your path to pharlap install dir"
+Quick smoke test that confirms PHaRLAP linked correctly and IRI is usable:
 
-6. export PYTHONPATH="Pylap_install_dir"
+```bash
+python3 -c "
+import importlib, math
+iri = importlib.import_module('pylap.iri2016')
+_, oarr = iri.iri2016(40.0, -90.0, 100.0, [2026, 1, 15, 18, 0], 100.0, 10.0, 50, {})
+foF2 = 8.98 * math.sqrt(max(oarr[0], 0)) / 1e6
+print(f'IRI OK — foF2={foF2:.1f} MHz, hmF2={oarr[1]:.0f} km')
+"
+```
 
-7. export LD_LIBRARY="/"YOUR PATH TO DIR"/l_comp_lib_2020.4.304_comp.for_redist/compilers_and_libraries_2020.4.304/linux/compiler/lib/intel64_lin" 
+Expected: `IRI OK — foF2=X.X MHz, hmF2=XXX km` (values depend on solar conditions at the queried date/time). If you see `ImportError` or a Fortran/link error, `PHARLAP_HOME` and `DIR_MODELS_REF_DAT` are probably unset or pointing at the wrong directory.
 
-8. export DIR_MODELS_REF_DAT="{PHARLAP_HOME}/dat"
+## IRI-2020 Ionosphere Model
 
-9. sudo apt-get install python3-tk python3-pil python3-pil.imagetk libqt5gui5 python3-pyqt5 
+PyLap now supports IRI-2020 via the `iri2020` Python package. This works **independently of PHaRLAP** for ionosphere generation.
 
-10. sudo apt-get install libxcb-randr0-dev libxcb-xtest0-dev libxcb-xinerama0-dev libxcb-shape0-dev libxcb-xkb-dev
+### Supported Profile Types
 
-11. source /home/{username}/bin/compilervars.sh intel64
+| Profile Type | Description |
+|--------------|-------------|
+| `'iri'` or `'iri2016'` | IRI-2016 model (default, requires PHaRLAP) |
+| `'iri2020'` | IRI-2020 model (standalone, no PHaRLAP needed) |
+| `'iri2012'` | IRI-2012 model (requires PHaRLAP) |
+| `'iri2007'` | IRI-2007 model (requires PHaRLAP) |
+| `'firi'` | IRI-2016 with FIRI D-region (requires PHaRLAP) |
 
-12. sudo apt-get install python3-pip
+### Example Usage
 
-13. pip3 install matplotlib numpy scipy qtpy
+```python
+from Ionosphere import gen_iono_grid_2d as gen_iono
 
-14. cd $PYTHONPATH 
+# Generate ionosphere using IRI-2020 (no PHaRLAP required)
+iono_pf_grid, iono_pf_grid_5, collision_freq, irreg, iono_te_grid = \
+    gen_iono.gen_iono_grid_2d(
+        origin_lat, origin_long, R12, UT, ray_bear,
+        max_range, num_range, range_inc, start_height,
+        height_inc, num_heights, kp, doppler_flag, 
+        'iri2020'  # Use IRI-2020
+    )
+```
 
-15. python3 setup.py install --user
+## Manual Setup
 
-16. Use Example folder files as templates to test the installation
+If you prefer not to use the setup script:
 
-    *Note* This is not a one time setup and will have to be redone if the terminal is closed out or if the code project is closed out *Note*
+```bash
+# 1. Set environment variables
+export PHARLAP_HOME="/path/to/pharlap_4.7.4"
+export PYTHONPATH="/path/to/PyLap"
+export DIR_MODELS_REF_DAT="${PHARLAP_HOME}/dat"
 
+# 2. Install system dependencies (Linux)
+sudo apt-get install python3-tk python3-pil python3-pil.imagetk \
+    libqt5gui5 python3-pyqt5 python3-venv gfortran
+sudo apt-get install libxcb-randr0-dev libxcb-xtest0-dev \
+    libxcb-xinerama0-dev libxcb-shape0-dev libxcb-xkb-dev
 
+# macOS equivalent
+# brew install gcc python-tk
 
-## for any questions or for help troubleshooting the install of PyLap please email Devin.diehl@scranton.edu ##
+# 3. Create virtual environment
+cd /path/to/PyLap
+python3 -m venv venv
+source venv/bin/activate
+
+# 4. Install Python packages
+pip install --upgrade pip
+pip install -r requirements.txt
+
+# 5. Build PyLap
+python3 setup.py install
+```
+
+## Re-running Setup
+
+To redo the setup:
+
+1. Edit `~/.bashrc` and remove PyLap environment variables
+2. Delete the venv: `rm -rf PyLap/venv`
+3. Run `setup.sh` again
+
+## Contact
+
+For questions or help troubleshooting, email: Devin.diehl@scranton.edu
diff --git a/modules/pylap/__init__.py b/modules/pylap/__init__.py
index 47b9900..c96deb9 100644
--- a/modules/pylap/__init__.py
+++ b/modules/pylap/__init__.py
@@ -1,16 +1,18 @@
-from pylap.abso_bg import abso_bg
-from pylap.dop_spread_eq import dop_spread_eq
-from pylap.ground_bs_loss import ground_bs_loss
-from pylap.ground_fs_loss import ground_fs_loss
-from pylap.igrf2007 import igrf2007
-from pylap.igrf2011 import igrf2011
-from pylap.igrf2016 import igrf2016
-from pylap.iri2007 import iri2007
-from pylap.iri2012 import iri2012
-from pylap.iri2016 import iri2016
-from pylap.irreg_strength import irreg_strength
-from pylap.nrlmsise00 import nrlmsise00
-from pylap.raytrace_2d import raytrace_2d
-from pylap.raytrace_2d_sp import raytrace_2d_sp
-from pylap.raytrace_3d import raytrace_3d
-from pylap.raytrace_3d_sp import raytrace_3d_sp
+# Soft imports — only expose modules that were actually compiled.
+# Modules requiring legacy IRI libs (iri2007, iri2012) are absent on PHaRLAP 4.7.4.
+_optional = [
+    'abso_bg', 'dop_spread_eq', 'ground_bs_loss', 'ground_fs_loss',
+    'igrf2007', 'igrf2011', 'igrf2016',
+    'iri2007', 'iri2012', 'iri2016',
+    'irreg_strength', 'nrlmsise00',
+    'raytrace_2d', 'raytrace_2d_sp',
+    'raytrace_3d', 'raytrace_3d_sp',
+]
+for _m in _optional:
+    try:
+        from importlib import import_module as _imp
+        _mod = _imp('pylap.' + _m)
+        globals()[_m] = getattr(_mod, _m)
+    except (ImportError, AttributeError):
+        pass
+del _optional, _m, _imp, _mod
diff --git a/modules/source/common/check_ref_data.c b/modules/source/common/check_ref_data.c
index cb029b7..2d0fc7e 100644
--- a/modules/source/common/check_ref_data.c
+++ b/modules/source/common/check_ref_data.c
@@ -8,6 +8,7 @@
 int check_iri2007(const char *);
 int check_iri2012(const char *);
 int check_iri2016(const char *);
+int check_iri2020(const char *);
 int check_land_sea(const char *);
 
 int check_file_exists(const char *);
@@ -24,6 +25,8 @@ int check_ref_data(const char *files_to_check)
     return check_iri2012(refdata_dir);
   } else if (strncmp(files_to_check, "iri2016", 7) == 0) {
     return check_iri2016(refdata_dir);
+  } else if (strncmp(files_to_check, "iri2020", 7) == 0) {
+    return check_iri2020(refdata_dir);
   } else if (strncmp(files_to_check, "land_sea", 8) == 0) {
     return check_land_sea(refdata_dir);
   }
@@ -134,6 +137,41 @@ int check_iri2016(const char *refdata_dir)
   return 1;
 }
 
+int check_iri2020(const char *refdata_dir)
+{
+  char filename[256];
+
+  snprintf(filename, 256, "%s/iri2020/igrf2025.dat", refdata_dir);
+  ASSERT_INT_NOMSG(check_file_exists(filename));
+
+  snprintf(filename, 256, "%s/iri2020/igrf2025s.dat", refdata_dir);
+  ASSERT_INT_NOMSG(check_file_exists(filename));
+
+  for (int year = 1945; year <= 2010; year += 5) {
+    snprintf(filename, 256, "%s/iri2020/dgrf%d.dat", refdata_dir, year);
+    ASSERT_INT_NOMSG(check_file_exists(filename));
+  }
+
+  snprintf(filename, 256, "%s/iri2020/apf107.dat", refdata_dir);
+  ASSERT_INT_NOMSG(check_file_exists(filename));
+
+  snprintf(filename, 256, "%s/iri2020/ig_rz.dat", refdata_dir);
+  ASSERT_INT_NOMSG(check_file_exists(filename));
+
+  for (int month = 1; month <= 12; month++) {
+    snprintf(filename, 256, "%s/iri2020/mcsat%d.dat", refdata_dir, month+10);
+    ASSERT_INT_NOMSG(check_file_exists(filename));
+
+    snprintf(filename, 256, "%s/iri2020/ccir%d.asc", refdata_dir, month+10);
+    ASSERT_INT_NOMSG(check_file_exists(filename));
+
+    snprintf(filename, 256, "%s/iri2020/ursi%d.asc", refdata_dir, month+10);
+    ASSERT_INT_NOMSG(check_file_exists(filename));
+  }
+
+  return 1;
+}
+
 int check_land_sea(const char *refdata_dir)
 {
   char filename[256];
diff --git a/modules/source/igrf2016.c b/modules/source/igrf2016.c
index 00e4ecc..4117100 100644
--- a/modules/source/igrf2016.c
+++ b/modules/source/igrf2016.c
@@ -1,6 +1,6 @@
 #include "../include/pharlap.h"
 
-extern void igrf2016_calc_(float *glat, float *glon, float *dec_year,
+extern void igrf2020_calc_(float *glat, float *glon, float *dec_year,
   float *height, float *dipole_moment, float *babs, float *bnorth, float *beast,
   float *bdown, float *dip, float *dec, float *dip_lat, float *l_value,
   int *l_value_code);
@@ -49,7 +49,7 @@ static PyObject *igrf2016(PyObject *self, PyObject *args)
   int l_value_code;
 
   /* Call subroutine. */
-  igrf2016_calc_(&latitude, &longitude, &dec_year, &height, &dipole_moment,
+  igrf2020_calc_(&latitude, &longitude, &dec_year, &height, &dipole_moment,
     &babs, &bnorth, &beast, &bdown, &dip, &dec, &dip_lat, &l_value,
     &l_value_code);
 
diff --git a/modules/source/iri2016.c b/modules/source/iri2016.c
index 07b89cf..ea34cbb 100644
--- a/modules/source/iri2016.c
+++ b/modules/source/iri2016.c
@@ -1,6 +1,6 @@
 #include "../include/pharlap.h"
 
-extern void iri2016_calc_(int *jf, int *jmag, float *glat, float *glon,
+extern void iri2020_calc_(int *jf, int *jmag, float *glat, float *glon,
   int *year, int *mmdd, float *dhour, float *heibeg, float *heiend,
   float *heistep, float *outf, float *oarr);
 
@@ -21,7 +21,7 @@ static PyObject *iri2016(PyObject *self, PyObject *args)
       &PyList_Type, &tm, &height_start, &height_step, &num_heights,
       &PyDict_Type, &iri_options));
 
-  ASSERT_NOMSG(check_ref_data("iri2016"));
+  ASSERT_NOMSG(check_ref_data("iri2020"));
 
   /* Parse UT. */
   int *ut = (int *)malloc(5 * sizeof(int));
@@ -562,7 +562,7 @@ if(obj!=NULL){
 
   int jmag = 0;
 
-  iri2016_calc_(jf, &jmag, &glat, &glon, &year, &mmdd, &dhour, &height_start,
+  iri2020_calc_(jf, &jmag, &glat, &glon, &year, &mmdd, &dhour, &height_start,
     &height_end, &height_step, outf, oarr);
 
   /*
diff --git a/modules/source/raytrace_2d.c b/modules/source/raytrace_2d.c
index d70b7cb..b0e3074 100644
--- a/modules/source/raytrace_2d.c
+++ b/modules/source/raytrace_2d.c
@@ -12,13 +12,13 @@ extern void raytrace_2d_(double *origin_lat, double *origin_lon, int *num_rays,
                          int *nhops_attempted, double *ray_state_vec_in, int *npts_in_ray, double *ray_state_vec_out,
                          double *elapsed_time);
 
-static void verifyIonoGrid(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
+static int verifyIonoGrid(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
                            PyArrayObject *irreg_grid);
 
 static PyObject *buildOutput(int num_rays, int *nhops_attempted, int *npts_in_ray, double *ray_data,
                              double *ray_path_data, double *freqs_data, double *elevs_data, int *ray_label, double *ray_state_vec_out);
 
-static void buidlIonoStruct(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
+static int buidlIonoStruct(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
                             PyArrayObject *irreg_grid, double height_start, double height_inc, double range_inc);
 
 static struct ionosphere_struct ionosphere;
@@ -30,14 +30,14 @@ static PyObject *raytrace_2d(PyObject *self, PyObject *args)
 
   Py_ssize_t num_args = PyTuple_Size(args);
 
-  ASSERT((num_args == 8 || num_args == 9 || num_args == 15 || num_args == 16),
+  ASSERT_INT((num_args == 8 || num_args == 9 || num_args == 15 || num_args == 16),
          PyExc_ValueError, "incorrect number of input arguments");
 
   int init_ionosphere = 1;
 
   if (num_args == 8 || num_args == 9)
   {
-    ASSERT(iono_exist_in_mem, PyExc_RuntimeError,
+    ASSERT_INT(iono_exist_in_mem, PyExc_RuntimeError,
            "the ionosphere has not been initialized");
 
     init_ionosphere = 0;
@@ -65,7 +65,7 @@ static PyObject *raytrace_2d(PyObject *self, PyObject *args)
   /* Ensure that `input_ray_state` is the correct size (if supplied). */
   if (num_args == 9 || num_args == 16)
   {
-    ASSERT((PyDict_Size(input_ray_state) == 9), PyExc_ValueError,
+    ASSERT_INT((PyDict_Size(input_ray_state) == 9), PyExc_ValueError,
            "incorrect number of fields in dictionary.");
 
     PyObject *items = PyDict_Values(input_ray_state);
@@ -74,41 +74,41 @@ static PyObject *raytrace_2d(PyObject *self, PyObject *args)
     {
       PyObject *item = PyList_GetItem(items, i);
 
-      ASSERT(PyArray_Check(item), PyExc_ValueError,
+      ASSERT_INT(PyArray_Check(item), PyExc_ValueError,
              "field value must be a numpy array");
-      ASSERT((PyArray_NDIM((PyArrayObject *)item) == 1), PyExc_ValueError,
+      ASSERT_INT((PyArray_NDIM((PyArrayObject *)item) == 1), PyExc_ValueError,
              "invalid shape for field");
 
       npy_intp *item_shape = PyArray_DIMS((PyArrayObject *)item);
 
-      ASSERT((item_shape[0] == elevs_shape[0]), PyExc_ValueError,
+      ASSERT_INT((item_shape[0] == elevs_shape[0]), PyExc_ValueError,
              "invalid size for field");
     }
   }
 
   /* Ensure that `elevs` and `freqs` are the same (and correct) size. */
-  ASSERT((PyArray_NDIM(elevs) == 1), PyExc_ValueError,
+  ASSERT_INT((PyArray_NDIM(elevs) == 1), PyExc_ValueError,
          "invalid shape for elevs");
-  ASSERT((PyArray_NDIM(freqs) == 1), PyExc_ValueError,
+  ASSERT_INT((PyArray_NDIM(freqs) == 1), PyExc_ValueError,
          "invalid shape for freqs");
 
   npy_intp *freqs_shape = PyArray_DIMS(elevs);
 
-  ASSERT((elevs_shape[0] == freqs_shape[0]), PyExc_ValueError,
+  ASSERT_INT((elevs_shape[0] == freqs_shape[0]), PyExc_ValueError,
          "shape of elevs and freqs must be identical");
 
   /* Ensure the ionosphere grids are valid. */
   if (init_ionosphere)
   {
-    verifyIonoGrid(iono_en_grid, iono_en_grid_5, collision_grid, irreg_grid);
+    if (!verifyIonoGrid(iono_en_grid, iono_en_grid_5, collision_grid, irreg_grid)) return NULL;
   }
 
   /* Ensure that `nhops` is valid (0 < nhops <= 50). */
-  ASSERT((nhops > 0 && nhops <= 50), PyExc_ValueError,
+  ASSERT_INT((nhops > 0 && nhops <= 50), PyExc_ValueError,
          "number of hops is invalid; must be within the range of 1 through 50");
 
   /* Ensure that `tol` is valid (can be an integer or list of 3 elements. */
-  ASSERT((
+  ASSERT_INT((
              (PyList_CheckExact(in_tol) && PyList_Size(in_tol) == 3) ||
              PyLong_CheckExact(in_tol) ||
              PyFloat_CheckExact(in_tol)),
@@ -140,8 +140,8 @@ static PyObject *raytrace_2d(PyObject *self, PyObject *args)
   /* Load ionosphere */
   if (init_ionosphere)
   {
-    buidlIonoStruct(iono_en_grid, iono_en_grid_5, collision_grid,
-                    irreg_grid, height_start, height_inc, range_inc);
+    if (!buidlIonoStruct(iono_en_grid, iono_en_grid_5, collision_grid,
+                    irreg_grid, height_start, height_inc, range_inc)) return NULL;
     iono_exist_in_mem = 1;
   }
 
@@ -165,14 +165,14 @@ static PyObject *raytrace_2d(PyObject *self, PyObject *args)
           "the field \"%s\" is missing from input_ray_state.",
           ray_state_fields[field]);
 
-      ASSERT((val != NULL), PyExc_ValueError, message);
+      ASSERT_INT((val != NULL), PyExc_ValueError, message);
 
       sprintf(
           message,
           "the field \"%s\" must be a NumPy array.",
           ray_state_fields[field]);
 
-      ASSERT(PyArray_Check(val), PyExc_ValueError, message);
+      ASSERT_INT(PyArray_Check(val), PyExc_ValueError, message);
 
       sprintf(
           message,
@@ -182,8 +182,8 @@ static PyObject *raytrace_2d(PyObject *self, PyObject *args)
       PyArrayObject *arr = (PyArrayObject *)val;
       arr = (PyArrayObject *)PyArray_Cast(arr, NPY_DOUBLE);
 
-      ASSERT((PyArray_NDIM(arr) == 1), PyExc_ValueError, message);
-      ASSERT((PyArray_DIMS(arr)[0] >= num_rays), PyExc_ValueError, message);
+      ASSERT_INT((PyArray_NDIM(arr) == 1), PyExc_ValueError, message);
+      ASSERT_INT((PyArray_DIMS(arr)[0] >= num_rays), PyExc_ValueError, message);
 
       for (npy_intp i = 0; i < num_rays; i++)
       {
@@ -206,7 +206,7 @@ static PyObject *raytrace_2d(PyObject *self, PyObject *args)
 
   int *nhops_attempted = (int *)malloc(num_rays * sizeof(int));
   int *npts_in_ray = (int *)malloc(num_rays * sizeof(int));
-  ray_data = (double *)malloc(19 * nhops * num_rays * sizeof(double));
+  ray_data = (double *)malloc(24 * nhops * num_rays * sizeof(double));
   ray_path_data = (double *)malloc(9 * MAX_POINTS_IN_RAY * num_rays * sizeof(double));
 
   int *ray_label = (int *)malloc(nhops * num_rays * sizeof(int));
@@ -240,17 +240,19 @@ static PyObject *raytrace_2d(PyObject *self, PyObject *args)
   return result;
 }
 
-static void verifyIonoGrid(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
+
+
+static int verifyIonoGrid(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
                            PyArrayObject *irreg_grid)
 {
 
-  ASSERT((PyArray_NDIM(iono_en_grid) == 2), PyExc_ValueError,
+  ASSERT_INT((PyArray_NDIM(iono_en_grid) == 2), PyExc_ValueError,
          "invalid shape for iono_en_grid");
-  ASSERT((PyArray_NDIM(iono_en_grid_5) == 2), PyExc_ValueError,
+  ASSERT_INT((PyArray_NDIM(iono_en_grid_5) == 2), PyExc_ValueError,
          "invalid shape for iono_en_grid_5");
-  ASSERT((PyArray_NDIM(collision_grid) == 2), PyExc_ValueError,
+  ASSERT_INT((PyArray_NDIM(collision_grid) == 2), PyExc_ValueError,
          "invalid shape for collision_grid");
-  ASSERT((PyArray_NDIM(irreg_grid) == 2), PyExc_ValueError,
+  ASSERT_INT((PyArray_NDIM(irreg_grid) == 2), PyExc_ValueError,
          "invalid shape for irreg_grid");
 
   npy_intp *iono_en_grid_shape = PyArray_DIMS(iono_en_grid);
@@ -258,32 +260,33 @@ static void verifyIonoGrid(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_g
   npy_intp *collision_grid_shape = PyArray_DIMS(collision_grid);
   npy_intp *irreg_grid_shape = PyArray_DIMS(irreg_grid);
 
-  ASSERT((
+  ASSERT_INT((
              iono_en_grid_shape[0] <= max_num_ht &&
              iono_en_grid_shape[1] <= max_num_rng),
          PyExc_ValueError, "iono_en_grid is too large");
-  ASSERT((
+  ASSERT_INT((
              iono_en_grid_5_shape[0] <= max_num_ht &&
              iono_en_grid_5_shape[1] <= max_num_rng),
          PyExc_ValueError, "iono_en_grid_5 is too large");
-  ASSERT((
+  ASSERT_INT((
              collision_grid_shape[0] <= max_num_ht &&
              collision_grid_shape[1] <= max_num_rng),
          PyExc_ValueError, "collision_grid is too large");
-  ASSERT((
+  ASSERT_INT((
              irreg_grid_shape[0] == 4 &&
              irreg_grid_shape[1] <= max_num_rng),
          PyExc_ValueError, "irreg_grid is not the correct size");
 
-  ASSERT((
+  ASSERT_INT((
              iono_en_grid_shape[0] == iono_en_grid_5_shape[0] &&
              iono_en_grid_shape[0] == collision_grid_shape[0]),
          PyExc_ValueError, "ionosphere grids have inconsistent row counts");
-  ASSERT((
+  ASSERT_INT((
              iono_en_grid_shape[1] == iono_en_grid_5_shape[1] &&
              iono_en_grid_shape[1] == collision_grid_shape[1] &&
              iono_en_grid_shape[1] == irreg_grid_shape[1]),
          PyExc_ValueError, "ionosphere grids have inconsistent column counts");
+  return 1;
 }
 
 static PyObject *buildOutput(int num_rays, int *nhops_attempted, int *npts_in_ray, double *ray_data,
@@ -330,7 +333,7 @@ static PyObject *buildOutput(int num_rays, int *nhops_attempted, int *npts_in_ra
       tmp = PyArray_ZEROS(1, nhops_dims, NPY_DOUBLE, 0);
       tmp_data = PyArray_DATA((PyArrayObject *)tmp);
 
-      stepmemcpyd(tmp_data, &ray_data[idx], num_rays * 9,
+      stepmemcpyd(tmp_data, &ray_data[idx], num_rays * 24,
                   nhops_attempted[ray_id]);
       PyDict_SetItemString(py_ray_data, rays_fields[field_id], tmp);
     }
@@ -390,7 +393,7 @@ static PyObject *buildOutput(int num_rays, int *nhops_attempted, int *npts_in_ra
   return PyTuple_Pack(3, py_rays, py_ray_paths, py_ray_states);
 }
 
-static void buidlIonoStruct(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
+static int buidlIonoStruct(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
                             PyArrayObject *irreg_grid, double height_start, double height_inc, double range_inc)
 {
 
@@ -439,8 +442,10 @@ static void buidlIonoStruct(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_
   ionosphere.HtMin = height_start;   // range_inc
   ionosphere.HtInc = height_inc;
   ionosphere.dRange = range_inc;
+  return 1;
 }
 
+
 static PyMethodDef methods[] = {
     {"raytrace_2d", raytrace_2d, METH_VARARGS, ""},
     {NULL, NULL, 0, NULL}};
diff --git a/modules/source/raytrace_3d.c b/modules/source/raytrace_3d.c
index e5c0c2e..fe95759 100644
--- a/modules/source/raytrace_3d.c
+++ b/modules/source/raytrace_3d.c
@@ -17,10 +17,10 @@ extern void raytrace_3d_(double *start_lat, double *start_long, double *start_he
 
 static PyObject *buildOutput(int num_rays, int *nhops_attempted, int *npts_in_ray, double *ray_data, double *bearings_data,
                              double *ray_path_data, double *freqs_data, double *elevs_data, int *ray_label, double *ray_state_vec_out);
-static void buildIonoStruct(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
+static int buildIonoStruct(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
                             PyArrayObject *Bx, PyArrayObject *By, PyArrayObject *Bz,
                             PyObject *iono_grid_parms, PyObject *geomag_grid_parms);
-static void verifyIono(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
+static int verifyIono(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
                        PyArrayObject *Bx, PyArrayObject *By, PyArrayObject *Bz,
                        PyObject *iono_grid_parms, PyObject *geomag_grid_parms);
 void stepmemcpyd(double *dest, double *src, int step, int num_vals);
@@ -43,6 +43,8 @@ static void clear_ionosphere(void)
   }
 }
 
+
+
 PyObject *raytrace_3d(PyObject *self, PyObject *args)
 {
   char message[256];
@@ -110,9 +112,9 @@ PyObject *raytrace_3d(PyObject *self, PyObject *args)
 
   if (init_ionosphere)
   {
-    verifyIono(iono_en_grid, iono_en_grid_5, collision_grid,
+    if (!verifyIono(iono_en_grid, iono_en_grid_5, collision_grid,
                Bx, By, Bz,
-               iono_grid_parms, geomag_grid_parms);
+               iono_grid_parms, geomag_grid_parms)) return NULL;;
   }
 
   /* Ensure that `start_lat` is valid (-90 < nhops <= 90). */
@@ -157,9 +159,9 @@ PyObject *raytrace_3d(PyObject *self, PyObject *args)
   if (init_ionosphere)
   {
 
-    buildIonoStruct(iono_en_grid, iono_en_grid_5, collision_grid,
+    if (!buildIonoStruct(iono_en_grid, iono_en_grid_5, collision_grid,
                     Bx, By, Bz,
-                    iono_grid_parms, geomag_grid_parms);
+                    iono_grid_parms, geomag_grid_parms)) return NULL;
     /* Now the ionosphere has been read in set the iono_exist_in_mem flag to
        indicate this for future raytrace calls */
     iono_exist_in_mem = 1;
@@ -400,93 +402,95 @@ static PyObject *buildOutput(int num_rays, int *nhops_attempted, int *npts_in_ra
   return PyTuple_Pack(3, py_rays, py_ray_paths, py_ray_states);
 }
 
-static void verifyIono(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
+static int verifyIono(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
                        PyArrayObject *Bx, PyArrayObject *By, PyArrayObject *Bz,
                        PyObject *iono_grid_parms, PyObject *geomag_grid_parms)
 {
-  ASSERT((PyArray_NDIM(iono_en_grid) == 3), PyExc_ValueError,
+  ASSERT_INT((PyArray_NDIM(iono_en_grid) == 3), PyExc_ValueError,
          "invalid shape for iono_en_grid");
-  ASSERT((PyArray_NDIM(iono_en_grid_5) == 3), PyExc_ValueError,
+  ASSERT_INT((PyArray_NDIM(iono_en_grid_5) == 3), PyExc_ValueError,
          "invalid shape for iono_en_grid_5");
-  ASSERT((PyArray_NDIM(collision_grid) == 3), PyExc_ValueError,
+  ASSERT_INT((PyArray_NDIM(collision_grid) == 3), PyExc_ValueError,
          "invalid shape for collision_grid");
 
   npy_intp *iono_en_grid_shape = PyArray_DIMS(iono_en_grid);
   npy_intp *iono_en_grid_5_shape = PyArray_DIMS(iono_en_grid_5);
   npy_intp *collision_grid_shape = PyArray_DIMS(collision_grid);
 
-  ASSERT((
+  ASSERT_INT((
              iono_en_grid_shape[0] <= max_num_lat &&
              iono_en_grid_shape[1] <= max_num_lon &&
              iono_en_grid_shape[2] <= max_num_ht),
          PyExc_ValueError, "iono_en_grid is too large");
-  ASSERT((
+  ASSERT_INT((
              iono_en_grid_5_shape[0] <= max_num_lat &&
              iono_en_grid_5_shape[1] <= max_num_lon &&
              iono_en_grid_5_shape[2] <= max_num_ht),
          PyExc_ValueError, "iono_en_grid_5 is too large");
-  ASSERT((
+  ASSERT_INT((
              collision_grid_shape[0] <= max_num_lat &&
              collision_grid_shape[1] <= max_num_lon &&
              collision_grid_shape[2] <= max_num_ht),
          PyExc_ValueError, "collision_grid is too large");
 
-  ASSERT((
+  ASSERT_INT((
              iono_en_grid_shape[0] == iono_en_grid_5_shape[0] &&
              iono_en_grid_shape[0] == collision_grid_shape[0]),
          PyExc_ValueError, "ionosphere grids have inconsistent row counts");
-  ASSERT((
+  ASSERT_INT((
              iono_en_grid_shape[1] == iono_en_grid_5_shape[1] &&
              iono_en_grid_shape[1] == collision_grid_shape[1]),
          PyExc_ValueError, "ionosphere grids have inconsistent column counts");
 
   // check B grid for shape and consistency
-  ASSERT((PyArray_NDIM(Bx) == 3), PyExc_ValueError,
+  ASSERT_INT((PyArray_NDIM(Bx) == 3), PyExc_ValueError,
          "invalid shape for Bx");
-  ASSERT((PyArray_NDIM(By) == 3), PyExc_ValueError,
+  ASSERT_INT((PyArray_NDIM(By) == 3), PyExc_ValueError,
          "invalid shape for By");
-  ASSERT((PyArray_NDIM(By) == 3), PyExc_ValueError,
+  ASSERT_INT((PyArray_NDIM(By) == 3), PyExc_ValueError,
          "invalid shape for Bz");
   npy_intp *Bx_shape = PyArray_DIMS(Bx);
   npy_intp *By_shape = PyArray_DIMS(By);
   npy_intp *Bz_shape = PyArray_DIMS(Bz);
 
-  ASSERT((
+  ASSERT_INT((
              Bx_shape[0] <= 101 &&
              Bx_shape[1] <= 101 &&
              Bx_shape[2] <= 201),
          PyExc_ValueError, "Bx is too large");
-  ASSERT((
+  ASSERT_INT((
              By_shape[0] <= 101 &&
              By_shape[1] <= 101 &&
              By_shape[2] <= 201),
          PyExc_ValueError, "By is too large");
-  ASSERT((
+  ASSERT_INT((
              Bz_shape[0] <= 101 &&
              Bz_shape[1] <= 101 &&
              Bz_shape[2] <= 201),
          PyExc_ValueError, "Bz is too large");
-  ASSERT((
+  ASSERT_INT((
              Bx_shape[0] == By_shape[0] &&
              By_shape[0] == Bz_shape[0]),
          PyExc_ValueError, "B grids have inconsistent row counts");
-  ASSERT((
+  ASSERT_INT((
              Bx_shape[1] == By_shape[1] &&
              By_shape[1] == Bz_shape[1]),
          PyExc_ValueError, "B grids have inconsistent column counts");
-  ASSERT((
+  ASSERT_INT((
              Bx_shape[2] == By_shape[2] &&
              By_shape[2] == Bz_shape[2]),
          PyExc_ValueError, "B grids have inconsistent height counts");
 
   // ensure the parms are valid
-  ASSERT((PyList_Size(iono_grid_parms) == 9), PyExc_ValueError,
+  ASSERT_INT((PyList_Size(iono_grid_parms) == 9), PyExc_ValueError,
          "invalid shape for iono_grid_parms");
-  ASSERT((PyList_Size(geomag_grid_parms) == 9), PyExc_ValueError,
+  ASSERT_INT((PyList_Size(geomag_grid_parms) == 9), PyExc_ValueError,
          "invalid shape for geomag_grid_parms");
+  return 1;
 }
 
-static void buildIonoStruct(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
+
+static int buildIonoStruct(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_grid_5, PyArrayObject *collision_grid,
                             PyArrayObject *Bx, PyArrayObject *By, PyArrayObject *Bz,
                             PyObject *iono_grid_parms, PyObject *geomag_grid_parms)
 {
@@ -523,24 +527,24 @@ static void buildIonoStruct(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_
   ptr_ionosphere->lat_min = PyFloat_AsDouble(PyList_GetItem(iono_grid_parms, 0)); //num of columns in ion_en-grid_5
   ptr_ionosphere->lat_inc = PyFloat_AsDouble(PyList_GetItem(iono_grid_parms, 1)); //num of rows in iono_en_grid_5
   // /* Ensure that `start_lat` is valid (-90 < nhops <= 90). */
-  ASSERT((ptr_ionosphere->lat_min >= -90.0 && ptr_ionosphere->lat_min <= 90.0), PyExc_ValueError,
+  ASSERT_INT((ptr_ionosphere->lat_min >= -90.0 && ptr_ionosphere->lat_min <= 90.0), PyExc_ValueError,
          "lat_min is invalid; must be within the range of -90 through 90");
 
   ptr_ionosphere->num_lat = PyFloat_AsDouble(PyList_GetItem(iono_grid_parms, 2)); //range_inc
   /* Ensure that `start_lon` is valid (-180 < nhops <= 180). */
-  ASSERT((grid_shape[0] == ptr_ionosphere->num_lat), PyExc_ValueError,
+  ASSERT_INT((grid_shape[0] == ptr_ionosphere->num_lat), PyExc_ValueError,
          "grid shape != num_lat");
 
   ptr_ionosphere->lat_max = ptr_ionosphere->lat_min +
                             (ptr_ionosphere->num_lat - 1) * ptr_ionosphere->lat_inc;
 
   ptr_ionosphere->lon_min = PyFloat_AsDouble(PyList_GetItem(iono_grid_parms, 3));
-  ASSERT((ptr_ionosphere->lon_min >= -180 && ptr_ionosphere->lon_min <= 180.0), PyExc_ValueError,
+  ASSERT_INT((ptr_ionosphere->lon_min >= -180 && ptr_ionosphere->lon_min <= 180.0), PyExc_ValueError,
          "lon_min must be within -180 and 180");
 
   ptr_ionosphere->lon_inc = PyFloat_AsDouble(PyList_GetItem(iono_grid_parms, 4));
   ptr_ionosphere->num_lon = PyFloat_AsDouble(PyList_GetItem(iono_grid_parms, 5));
-  ASSERT((grid_shape[1] == ptr_ionosphere->num_lon), PyExc_ValueError,
+  ASSERT_INT((grid_shape[1] == ptr_ionosphere->num_lon), PyExc_ValueError,
          "grid shape[1] != num_lon");
 
   ptr_ionosphere->lon_max = ptr_ionosphere->lon_min +
@@ -548,7 +552,7 @@ static void buildIonoStruct(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_
   ptr_ionosphere->ht_min = PyFloat_AsDouble(PyList_GetItem(iono_grid_parms, 6));
   ptr_ionosphere->ht_inc = PyFloat_AsDouble(PyList_GetItem(iono_grid_parms, 7));
   ptr_ionosphere->num_ht = PyFloat_AsDouble(PyList_GetItem(iono_grid_parms, 8));
-  ASSERT((grid_shape[2] == ptr_ionosphere->num_ht), PyExc_ValueError,
+  ASSERT_INT((grid_shape[2] == ptr_ionosphere->num_ht), PyExc_ValueError,
          "grid shape[2] != num_ht");
 
   ptr_ionosphere->ht_max = ptr_ionosphere->ht_min +
@@ -570,24 +574,24 @@ static void buildIonoStruct(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_
     }
   }
   geomag_field.lat_min = PyFloat_AsDouble(PyList_GetItem(geomag_grid_parms, 0));
-  ASSERT((geomag_field.lat_min >= -90.0 && geomag_field.lat_min <= 90), PyExc_ValueError,
+  ASSERT_INT((geomag_field.lat_min >= -90.0 && geomag_field.lat_min <= 90), PyExc_ValueError,
          "The start latitude of the input geomagnetic field grid must be in the range -90 to 90 degrees.");
 
   geomag_field.lat_inc = PyFloat_AsDouble(PyList_GetItem(geomag_grid_parms, 1));
   geomag_field.num_lat = PyFloat_AsDouble(PyList_GetItem(geomag_grid_parms, 2));
-  ASSERT((dims[0] == geomag_field.num_lat), PyExc_ValueError,
+  ASSERT_INT((dims[0] == geomag_field.num_lat), PyExc_ValueError,
          "The number of latitudes in the input geomagnetic field grid does not match the input grid parameter.");
 
   geomag_field.lat_max = geomag_field.lat_min +
                          (geomag_field.num_lat - 1) * geomag_field.lat_inc;
   geomag_field.lon_min = PyFloat_AsDouble(PyList_GetItem(geomag_grid_parms, 3));
 
-  ASSERT((geomag_field.lon_min >= -180.0 && geomag_field.lon_min <= 180.0), PyExc_ValueError,
+  ASSERT_INT((geomag_field.lon_min >= -180.0 && geomag_field.lon_min <= 180.0), PyExc_ValueError,
          "he start longitude of the input geomagnetic field grid must be in the range -180 to 180 degrees.");
 
   geomag_field.lon_inc = PyFloat_AsDouble(PyList_GetItem(geomag_grid_parms, 4));
   geomag_field.num_lon = PyFloat_AsDouble(PyList_GetItem(geomag_grid_parms, 5));
-  ASSERT((dims[1] == geomag_field.num_lon), PyExc_ValueError,
+  ASSERT_INT((dims[1] == geomag_field.num_lon), PyExc_ValueError,
          "The number of longitudes in the input geomagnetic field grid does not match the input grid parameter.");
 
   geomag_field.lon_max = geomag_field.lon_min +
@@ -595,12 +599,14 @@ static void buildIonoStruct(PyArrayObject *iono_en_grid, PyArrayObject *iono_en_
   geomag_field.ht_min = PyFloat_AsDouble(PyList_GetItem(geomag_grid_parms, 6));
   geomag_field.ht_inc = PyFloat_AsDouble(PyList_GetItem(geomag_grid_parms, 7));
   geomag_field.num_ht = PyFloat_AsDouble(PyList_GetItem(geomag_grid_parms, 8));
-  ASSERT((dims[2] == geomag_field.num_ht), PyExc_ValueError,
+  ASSERT_INT((dims[2] == geomag_field.num_ht), PyExc_ValueError,
          "The number of heights in the input geomagnetic field grid does not match the input grid parameter.");
 
   geomag_field.ht_max = geomag_field.ht_min +
                         (geomag_field.num_ht - 1) * geomag_field.ht_inc;
+  return 1;
 }
+
 static PyMethodDef methods[] = {
     {"raytrace_3d", raytrace_3d, METH_VARARGS, ""},
     {NULL, NULL, 0, NULL}};
diff --git a/requirements.txt b/requirements.txt
index 56735a8..207c245 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,12 +1,14 @@
+# Python packages (install with pip)
 scipy
 matplotlib
-python3-tk
-pthon3-phil 
-python3-phil.imagetk
+numpy
 qtpy
 PyQt5
-libxcb-randr0-dev 
-libxcb-xtest0-dev 
-libxcb-xinerama0-dev 
-libxcb-shape0-dev 
-libxcb-xkb-dev
\ No newline at end of file
+xarray
+
+# IRI-2020 must be installed from GitHub (requires gfortran):
+# pip install git+https://github.com/space-physics/iri2020.git
+
+# System packages (install with apt-get, not pip):
+# sudo apt-get install python3-tk python3-pil python3-pil.imagetk gfortran
+# sudo apt-get install libxcb-randr0-dev libxcb-xtest0-dev libxcb-xinerama0-dev libxcb-shape0-dev libxcb-xkb-dev
\ No newline at end of file
diff --git a/setup.py b/setup.py
index af4d300..524ab77 100644
--- a/setup.py
+++ b/setup.py
@@ -1,87 +1,142 @@
-    
 import glob
 import os
 import platform
+import struct
 
-#import setuptools
 import numpy as np
-from distutils.core import setup
-from distutils.core import Extension
-
-
-if platform.system() not in ['Linux']:
-    raise OSError('Unrecognized or unsupported operating system.')
-
-#
-# Configure PHaRLAP.
-#
+try:
+    from setuptools import setup, Extension
+except ImportError:
+    from distutils.core import setup
+    from distutils.core import Extension
+
+
+# ---------------------------------------------------------------------------
+# Platform detection
+# ---------------------------------------------------------------------------
+_sys = platform.system()
+_machine = platform.machine()
+
+# PHaRLAP 4.7.4 ships GCC-compiled static libraries on all platforms.
+# We need gfortran runtime on both Linux and macOS.
+import subprocess, shutil
+
+if _sys == 'Linux':
+    _lib_subdir = 'linux'
+elif _sys == 'Darwin':
+    _lib_subdir = 'maca' if _machine == 'arm64' else 'maci'
+else:
+    raise OSError(f'Unsupported operating system: {_sys}')
+
+# Locate gfortran runtime library directory.
+_extra_lib_dirs = []
+_gfc = shutil.which('gfortran')
+if _gfc:
+    _libname = 'libgfortran.dylib' if _sys == 'Darwin' else 'libgfortran.so'
+    _gfc_libdir = subprocess.check_output(
+        [_gfc, f'-print-file-name={_libname}'],
+        text=True
+    ).strip()
+    _gfc_libdir = os.path.dirname(_gfc_libdir)
+    if _gfc_libdir:
+        _extra_lib_dirs = [_gfc_libdir]
+elif _sys == 'Darwin':
+    _extra_lib_dirs = ['/opt/homebrew/lib/gcc/current']
+
+_extra_libs = ['gfortran', 'gomp']
+
+# ---------------------------------------------------------------------------
+# Locate PHaRLAP
+# ---------------------------------------------------------------------------
 if 'PHARLAP_HOME' not in os.environ:
-    raise OSError('The environment variable "PHARLAP_HOME" must be defined.')
+    raise OSError('Set PHARLAP_HOME to the PHaRLAP install directory before building.')
 
-if 'LD_LIBRARY' not in os.environ:
-    raise OSError('The environment variable "LD_LIBRARY" must be defined.')
-
-if 'PYTHONPATH' not in os.environ:
-    raise OSError('The environment variable "PYTHONPATH" must be defined.')
-
-pharlap_path = os.getenv('PHARLAP_HOME')
-intel_path  = os.getenv('LD_LIBRARY')
-py_path = os.getenv('PYTHONPATH')
+pharlap_path = os.environ['PHARLAP_HOME']
+if not os.path.isdir(pharlap_path):
+    raise OSError(f'PHARLAP_HOME does not exist: {pharlap_path}')
 
-if not os.path.isdir(py_path):
-    raise OSError('The environment variable "PYTHONPATH" is invalid.')
+pharlap_include_path = os.path.join(pharlap_path, 'src', 'C')
+pharlap_lib_path = os.path.join(pharlap_path, 'lib', _lib_subdir)
 
-if not os.path.isdir(pharlap_path):
-    raise OSError('The environment variable "PHARLAP_HOME" is invalid.')
+if not os.path.isdir(pharlap_lib_path):
+    raise OSError(f'PHaRLAP library path not found: {pharlap_lib_path}')
 
-if not os.path.isdir(intel_path):
-    raise OSError('The environment variable "LD_LIBRARY" is invalid.')
+# No Intel Fortran needed — PHaRLAP 4.7.4 libs are GCC-compiled on all platforms.
 
-pharlap_include_path = os.path.join(pharlap_path, 'src', 'C')
-pharlap_lib_path = os.path.join(pharlap_path, 'lib', 'linux')
-# Define native modules.
-#
+# PHaRLAP 4.7.4 ships libiri2020 (consolidated; replaces iri2007/2012/2016).
+# Modules that needed legacy IRI versions are skipped when those libs are absent.
+_available_libs = {
+    os.path.splitext(f)[0].replace('lib', '')
+    for f in os.listdir(pharlap_lib_path)
+    if f.startswith('lib')
+}
 
+# ---------------------------------------------------------------------------
+# Module definitions
+# ---------------------------------------------------------------------------
 native_modules = []
 COMMON_GLOB = glob.glob('modules/source/common/*.c')
 
+
+def _libs(*base):
+    """Return the base PHaRLAP libs plus platform math lib."""
+    return list(base) + ['m']
+
+
 def create_module(name, libraries):
+    """Register a C-extension module if all required native libs are present."""
+    missing = [lib for lib in libraries if lib not in _available_libs
+               and lib not in ('m',)]
+    if missing:
+        print(f'  SKIP  pylap.{name}  (missing PHaRLAP libs: {missing})')
+        return
+    src = os.path.join('modules', 'source', name + '.c')
+    if not os.path.exists(src):
+        print(f'  SKIP  pylap.{name}  (source not found: {src})')
+        return
     native_modules.append(Extension(
         'pylap.' + name,
-        sources=['modules/source/' + name + '.c'] + COMMON_GLOB,
-        include_dirs=[np.get_include(), pharlap_include_path, "include"],
-        library_dirs=[pharlap_lib_path, intel_path],
-        libraries=libraries))
-    pass
-
-create_module('abso_bg', ['propagation', 'maths', 'iri2016', 'ifcore', 'imf', 'irc', 'svml'])
-create_module('dop_spread_eq', ['propagation', 'iri2016', 'ifcore', 'imf', 'irc', 'svml'])
-create_module('ground_bs_loss', ['propagation', 'ifcore', 'imf'])
-create_module('ground_fs_loss', ['propagation', 'ifcore', 'imf'])
-create_module('igrf2007', ['maths', 'iri2007', 'ifcore', 'imf', 'irc'])
-create_module('igrf2011', ['maths', 'iri2012', 'ifcore', 'imf', 'irc'])
-create_module('igrf2016', ['maths', 'iri2016', 'ifcore', 'imf', 'irc', 'svml'])
-create_module('iri2007', ['iri2007', 'ifcore', 'imf', 'irc', 'svml'])
-create_module('iri2012', ['iri2012', 'ifcore', 'imf', 'irc', 'svml'])
-create_module('iri2016', ['iri2016', 'ifcore', 'imf', 'irc', 'svml'])
-create_module('irreg_strength', ['propagation', 'iri2016', 'ifcore', 'imf', 'irc', 'svml'])
-create_module('nrlmsise00', ['maths', 'iri2016', 'ifcore', 'imf', 'irc', 'svml'])
-create_module('raytrace_2d', ['propagation', 'maths', 'ifcore','imf','iomp5'])
-create_module('raytrace_2d_sp', ['propagation', 'maths', 'ifcore', 'imf', 'iomp5'])
-create_module('raytrace_3d', ['propagation', 'maths', 'ifcore', 'imf', 'iomp5'])
-create_module('raytrace_3d_sp', ['propagation', 'maths', 'ifcore', 'imf', 'iomp5'])
-
-#
-# Create module.
-#
-#with open('requirements.txt') as f:
- #   requirements = f.read().splitlines()
+        sources=[src] + COMMON_GLOB,
+        include_dirs=[np.get_include(), pharlap_include_path,
+                      os.path.join('modules', 'include')],
+        library_dirs=[pharlap_lib_path] + _extra_lib_dirs,
+        libraries=libraries + _extra_libs,
+        extra_compile_args=["-Wno-error=incompatible-pointer-types"],
+    ))
+    print(f'  BUILD pylap.{name}')
+
+
+# PHaRLAP 4.7.4: IRI consolidated into iri2020.  Modules that required
+# the legacy iri2007/iri2012/iri2016 archives are updated accordingly;
+# those that only existed as thin wrappers for the legacy models are skipped.
+print(f'\nConfiguring pyLAP for PHaRLAP on {_sys}/{_machine} (lib/{_lib_subdir})\n')
+
+create_module('abso_bg',       _libs('propagation', 'maths', 'iri2020'))
+create_module('dop_spread_eq', _libs('propagation', 'iri2020'))
+create_module('ground_bs_loss',_libs('propagation'))
+create_module('ground_fs_loss',_libs('propagation'))
+# igrf2007/igrf2011 need legacy iri2007/iri2012 — skipped on PHaRLAP 4.7.4
+create_module('igrf2007',      _libs('maths', 'iri2007'))
+create_module('igrf2011',      _libs('maths', 'iri2012'))
+create_module('igrf2016',      _libs('maths', 'iri2020'))
+# Standalone IRI wrappers — skipped for legacy versions not in 4.7.4
+create_module('iri2007',       _libs('iri2007'))
+create_module('iri2012',       _libs('iri2012'))
+create_module('iri2016',       _libs('iri2020'))
+create_module('irreg_strength',_libs('propagation', 'iri2020'))
+create_module('nrlmsise00',    _libs('maths', 'iri2020'))
+create_module('raytrace_2d',   _libs('propagation', 'maths'))
+create_module('raytrace_2d_sp',_libs('propagation', 'maths'))
+create_module('raytrace_3d',   _libs('propagation', 'maths'))
+create_module('raytrace_3d_sp',_libs('propagation', 'maths'))
+
+print()
 
 setup(
-        name='pylap',
-        packages=['modules/pylap'],
-        #install_requires=requirements,
-        version='0.1.0-alpha',
-        description='A numpy-compatible Python 3 wrapper for the PHaRLAP ionospheric raytracer',
-        ext_modules=native_modules
+    name='pylap',
+    version='0.1.0-alpha',
+    description='A numpy-compatible Python 3 wrapper for the PHaRLAP ionospheric raytracer',
+    packages=['pylap'],
+    package_dir={'pylap': 'modules/pylap'},
+    ext_modules=native_modules,
 )
diff --git a/setup.sh b/setup.sh
index aa356ac..602bd89 100644
--- a/setup.sh
+++ b/setup.sh
@@ -1,17 +1,32 @@
+#!/bin/bash
+# PyLap Setup Script
+# Uses a Python virtual environment for package management
 
-# check if the program is already setup
+set -e
 
-if [[ ${PHARLAP_HOME} = "" ]]; then
+SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+VENV_DIR="${SCRIPT_DIR}/venv"
 
+# check if the program is already setup
+if [[ ${PHARLAP_HOME} = "" ]]; then
 
 echo "First Time Setup"
-echo "Please enter the filepath of the project"
-echo "Example - /home/devindiehl/Pylap_project"
+echo "Please enter the filepath of the project directory"
+echo "Example - /home/username/pylap_project"
 read x
 echo ${x}
 
-export PHARLAP_HOME="${x}/pharlap_4.5.1"
-echo "export PHARLAP_HOME=${x}/pharlap_4.5.1" >> ~/.bashrc
+# Detect PHaRLAP version (look for pharlap_* directory)
+PHARLAP_DIR=$(find "${x}" -maxdepth 1 -type d -name "pharlap_*" 2>/dev/null | head -1)
+if [[ -n "${PHARLAP_DIR}" ]]; then
+    export PHARLAP_HOME="${PHARLAP_DIR}"
+    echo "Found PHaRLAP at: ${PHARLAP_HOME}"
+else
+    echo "Warning: PHaRLAP directory not found in ${x}"
+    echo "Please ensure PHaRLAP is installed before running raytracing examples"
+    export PHARLAP_HOME="${x}/pharlap_4.7.4"
+fi
+echo "export PHARLAP_HOME=${PHARLAP_HOME}" >> ~/.bashrc
 
 if [[ -d "${x}/PyLap" ]] ;
 then
@@ -27,37 +42,49 @@ else
 echo "Pylap filepath not found"
 fi
 
-export LD_LIBRARY="${x}/l_comp_lib_2020.4.304_comp.for_redist/compilers_and_libraries_2020.4.304/linux/compiler/lib/intel64_lin"
-echo "export LD_LIBRARY=${x}/l_comp_lib_2020.4.304_comp.for_redist/compilers_and_libraries_2020.4.304/linux/compiler/lib/intel64_lin" >> ~/.bashrc
-
-export DIR_MODELS_REF_DAT="${x}/pharlap_4.5.1/dat"
-echo "export DIR_MODELS_REF_DAT=${x}/pharlap_4.5.1/dat">> ~/.bashrc
-
-source /home/$USER/bin/compilervars.sh intel64
-echo "source /home/$USER/bin/compilervars.sh intel64" >>  ~/.bashrc
-
-
-# echo "python3 setup.py install --user" >> ~/.bashrc
-
-# export PHARLAP_HOME="/home/devindiehl/Pylap_project/PHARLAP/pharlap_4.5.1"
-
-# export PYTHONPATH="/home/devindiehl/Pylap_project/PyLap"
-
-# export LD_LIBRARY="/home/devindiehl/Pylap_project/l_comp_lib_2020.4.304_comp.for_redist/compilers_and_libraries_2020.4.304/linux/compiler/lib/intel64_lin" 
-
-# export DIR_MODELS_REF_DAT="/home/devindiehl/Pylap_project/PHARLAP/pharlap_4.5.1/dat"
-
-
-sudo apt-get install python3-tk python3-pil python3-pil.imagetk libqt5gui5 python3-pyqt5 
-
-sudo apt-get install libxcb-randr0-dev libxcb-xtest0-dev libxcb-xinerama0-dev libxcb-shape0-dev libxcb-xkb-dev
-
-sudo apt-get install python3-pip
-
-pip3 install matplotlib numpy scipy qtpy
+export DIR_MODELS_REF_DAT="${PHARLAP_HOME}/dat"
+echo "export DIR_MODELS_REF_DAT=${DIR_MODELS_REF_DAT}" >> ~/.bashrc
+
+# Install system dependencies (gfortran provides the Fortran runtime — Intel Fortran is no longer required)
+echo "Installing system dependencies..."
+sudo apt-get update
+sudo apt-get install -y python3-tk python3-pil python3-pil.imagetk libqt5gui5 python3-pyqt5 
+sudo apt-get install -y libxcb-randr0-dev libxcb-xtest0-dev libxcb-xinerama0-dev libxcb-shape0-dev libxcb-xkb-dev
+sudo apt-get install -y python3-pip python3-venv gfortran
+
+# Create and activate virtual environment
+echo "Creating Python virtual environment at ${VENV_DIR}..."
+python3 -m venv "${VENV_DIR}"
+source "${VENV_DIR}/bin/activate"
+
+# Install Python dependencies in venv
+echo "Installing Python packages in virtual environment..."
+pip install --upgrade pip
+pip install -r "${SCRIPT_DIR}/requirements.txt"
+
+# Install IRI-2020 from GitHub (requires gfortran)
+echo "Installing IRI-2020 from GitHub..."
+pip install git+https://github.com/space-physics/iri2020.git
+
+# Build and install PyLap
+echo "Building PyLap..."
+python3 setup.py install
+
+# Add venv activation to bashrc
+echo "# PyLap virtual environment" >> ~/.bashrc
+echo "alias pylap-activate='source ${VENV_DIR}/bin/activate'" >> ~/.bashrc
+echo ""
+echo "Setup complete!"
+echo "To activate the PyLap environment, run: source ${VENV_DIR}/bin/activate"
+echo "Or use the alias: pylap-activate"
 
-python3 setup.py install --user
 else
 echo "Pylap is already setup"
-echo "if you wish to redo the setup enter the command 'nano .bashrc' in the home directory, delete the filepaths and then run the setup.sh script again"
+echo "To activate the virtual environment: source ${SCRIPT_DIR}/venv/bin/activate"
+echo "Or use the alias: pylap-activate"
+echo ""
+echo "If you wish to redo the setup:"
+echo "1. Edit ~/.bashrc and remove the PyLap environment variables"
+echo "2. Delete the venv directory: rm -rf ${SCRIPT_DIR}/venv"
+echo "3. Run this script again"
 fi