change stop to mesa_error, real to real(dp) in colors

Author: pmocz

colors/private/colors_utils.f90

@@ -20,6 +20,7 @@
 module colors_utils
    use const_def, only: dp, strlen, mesa_dir
    use colors_def, only: Colors_General_Info
+   use utils_lib, only: mesa_error
 
    implicit none
 
@@ -39,7 +40,7 @@ subroutine dilute_flux(surface_flux, R, d, calibrated_flux)
       ! Check that the output array has the same size as the input
       if (size(calibrated_flux) /= size(surface_flux)) then
          print *, "Error in dilute_flux: Output array must have the same size as input array."
-         stop 1
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       ! Apply the dilution factor (R/d)^2 to each element
@@ -59,25 +60,25 @@ subroutine trapezoidal_integration(x, y, result)
       real(dp), intent(out) :: result
 
       integer :: i, n
-      REAL :: sum
+      real(dp) :: sum
 
       n = size(x)
-      sum = 0.0
+      sum = 0.0_dp
 
       ! Validate input sizes
       if (size(x) /= size(y)) then
          print *, "Error: x and y arrays must have the same size."
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       if (size(x) < 2) then
          print *, "Error: x and y arrays must have at least 2 elements."
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       ! Perform trapezoidal integration
       do i = 1, n - 1
-         sum = sum + 0.5*(x(i + 1) - x(i))*(y(i + 1) + y(i))
+         sum = sum + 0.5_dp*(x(i + 1) - x(i))*(y(i + 1) + y(i))
       end do
 
       result = sum
@@ -96,15 +97,15 @@ subroutine simpson_integration(x, y, result)
       ! Validate input sizes
       if (size(x) /= size(y)) then
          print *, "Error: x and y arrays must have the same size."
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       if (size(x) < 2) then
          print *, "Error: x and y arrays must have at least 2 elements."
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
-      ! Perform adaptive Simpson’s rule
+      ! Perform adaptive Simpson's rule
       do i = 1, n - 2, 2
          h1 = x(i + 1) - x(i)       ! Step size for first interval
          h2 = x(i + 2) - x(i + 1)     ! Step size for second interval
@@ -139,12 +140,12 @@ subroutine romberg_integration(x, y, result)
       ! Validate input sizes
       if (size(x) /= size(y)) then
          print *, "Error: x and y arrays must have the same size."
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       if (n < 2) then
          print *, "Error: x and y arrays must have at least 2 elements."
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       allocate (R(m))
@@ -192,14 +193,14 @@ subroutine load_vega_sed(filepath, wavelengths, flux)
       open (unit, file=trim(filepath), status='OLD', action='READ', iostat=status)
       if (status /= 0) then
          print *, "Error: Could not open Vega SED file ", trim(filepath)
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       ! Skip header line
       read (unit, '(A)', iostat=status) line
       if (status /= 0) then
          print *, "Error: Could not read header from Vega SED file ", trim(filepath)
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       ! Count the number of data lines
@@ -237,21 +238,21 @@ subroutine load_filter(directory, filter_wavelengths, filter_trans)
 
       character(len=512) :: line
       integer :: unit, n_rows, status, i
-      REAL :: temp_wavelength, temp_trans
+      real(dp) :: temp_wavelength, temp_trans
 
       ! Open the file
       unit = 20
       open (unit, file=trim(directory), status='OLD', action='READ', iostat=status)
       if (status /= 0) then
          print *, "Error: Could not open file ", trim(directory)
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       ! Skip header line
       read (unit, '(A)', iostat=status) line
       if (status /= 0) then
          print *, "Error: Could not read the file", trim(directory)
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       ! Count rows in the file
@@ -272,7 +273,7 @@ subroutine load_filter(directory, filter_wavelengths, filter_trans)
          read (unit, '(A)', iostat=status) line
          if (status /= 0) then
             print *, "Error: Could not rewind file", trim(directory)
-            stop
+            call mesa_error(__FILE__, __LINE__)
          end if
          if (line(1:1) /= "#") exit
       end do
@@ -313,14 +314,14 @@ subroutine load_lookup_table(lookup_file, lookup_table, out_file_names, &
       open (unit, file=lookup_file, status='old', action='read', iostat=status)
       if (status /= 0) then
          print *, "Error: Could not open file", lookup_file
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       ! Read header line
       read (unit, '(A)', iostat=status) line
       if (status /= 0) then
          print *, "Error: Could not read header line"
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       call split_line(line, delimiter, headers)
@@ -474,15 +475,15 @@ subroutine load_sed(directory, index, wavelengths, flux)
       open (unit, file=trim(directory), status='OLD', action='READ', iostat=status)
       if (status /= 0) then
          print *, "Error: Could not open file ", trim(directory)
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       ! Skip header lines
       do
          read (unit, '(A)', iostat=status) line
          if (status /= 0) then
             print *, "Error: Could not read the file", trim(directory)
-            stop
+            call mesa_error(__FILE__, __LINE__)
          end if
          if (line(1:1) /= "#") exit
       end do
@@ -505,7 +506,7 @@ subroutine load_sed(directory, index, wavelengths, flux)
          read (unit, '(A)', iostat=status) line
          if (status /= 0) then
             print *, "Error: Could not rewind file", trim(directory)
-            stop
+            call mesa_error(__FILE__, __LINE__)
          end if
          if (line(1:1) /= "#") exit
       end do
@@ -584,7 +585,7 @@ subroutine read_strings_from_file(colors_settings, strings, n, ierr)
       if (status /= 0) then
          ierr = -1
          print *, "Error: Could not open file", filename
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       do

colors/private/knn_interp.f90

@@ -24,6 +24,7 @@
 module knn_interp
    use const_def, only: dp
    use colors_utils, only: dilute_flux, load_sed
+   use utils_lib, only: mesa_error
    implicit none
 
    private
@@ -80,8 +81,8 @@ subroutine construct_sed_knn(teff, log_g, metallicity, R, d, file_names, &
          distances(i) = sqrt((lu_teff(closest_indices(i)) - teff)**2 + &
                              (lu_logg(closest_indices(i)) - log_g)**2 + &
                              (lu_meta(closest_indices(i)) - metallicity)**2)
-         if (distances(i) == 0.0) distances(i) = 1.0d-10  ! Prevent division by zero
-         weights(i) = 1.0/distances(i)
+         if (distances(i) == 0.0_dp) distances(i) = 1.0d-10  ! Prevent division by zero
+         weights(i) = 1.0_dp/distances(i)
       end do
 
       ! Normalize weights
@@ -91,7 +92,7 @@ subroutine construct_sed_knn(teff, log_g, metallicity, R, d, file_names, &
       ! Allocate output arrays
       allocate (wavelengths(n_points), fluxes(n_points))
       wavelengths = common_wavelengths
-      fluxes = 0.0
+      fluxes = 0.0_dp
 
       ! Perform weighted combination of the model fluxes (still at the stellar surface)
       do i = 1, 4
@@ -139,15 +140,15 @@ subroutine get_closest_stellar_models(teff, log_g, metallicity, lu_teff, &
       ! Allocate and scale lookup table values
       allocate (scaled_lu_teff(n), scaled_lu_logg(n), scaled_lu_meta(n))
 
-      if (teff_max - teff_min > 0.00) then
+      if (teff_max - teff_min > 0.0_dp) then
          scaled_lu_teff = (lu_teff - teff_min)/(teff_max - teff_min)
       end if
 
-      if (logg_max - logg_min > 0.00) then
+      if (logg_max - logg_min > 0.0_dp) then
          scaled_lu_logg = (lu_logg - logg_min)/(logg_max - logg_min)
       end if
 
-      if (meta_max - meta_min > 0.00) then
+      if (meta_max - meta_min > 0.0_dp) then
          scaled_lu_meta = (lu_meta - meta_min)/(meta_max - meta_min)
       end if
 
@@ -156,38 +157,35 @@ subroutine get_closest_stellar_models(teff, log_g, metallicity, lu_teff, &
       norm_logg = (log_g - logg_min)/(logg_max - logg_min)
       norm_meta = (metallicity - meta_min)/(meta_max - meta_min)
 
+      ! Detect dummy axes once (outside the loop)
+      use_teff_dim = .not. (all(lu_teff == 0.0_dp) .or. all(lu_teff == 999.0_dp) .or. all(lu_teff == -999.0_dp))
+      use_logg_dim = .not. (all(lu_logg == 0.0_dp) .or. all(lu_logg == 999.0_dp) .or. all(lu_logg == -999.0_dp))
+      use_meta_dim = .not. (all(lu_meta == 0.0_dp) .or. all(lu_meta == 999.0_dp) .or. all(lu_meta == -999.0_dp))
+
       ! Find closest models
       do i = 1, n
-         teff_dist = 0.0
-         logg_dist = 0.0
-         meta_dist = 0.0
+         teff_dist = 0.0_dp
+         logg_dist = 0.0_dp
+         meta_dist = 0.0_dp
 
-         if (teff_max - teff_min > 0.00) then
+         if (teff_max - teff_min > 0.0_dp) then
             teff_dist = scaled_lu_teff(i) - norm_teff
          end if
 
-         if (logg_max - logg_min > 0.00) then
+         if (logg_max - logg_min > 0.0_dp) then
             logg_dist = scaled_lu_logg(i) - norm_logg
          end if
 
-         if (meta_max - meta_min > 0.00) then
+         if (meta_max - meta_min > 0.0_dp) then
             meta_dist = scaled_lu_meta(i) - norm_meta
          end if
 
-         ! Detect dummy axes once
-
-         use_teff_dim = .not. (all(lu_teff == 0.0_dp) .or. all(lu_teff == 999.0_dp) .or. all(lu_teff == -999.0_dp))
-         use_logg_dim = .not. (all(lu_logg == 0.0_dp) .or. all(lu_logg == 999.0_dp) .or. all(lu_logg == -999.0_dp))
-         use_meta_dim = .not. (all(lu_meta == 0.0_dp) .or. all(lu_meta == 999.0_dp) .or. all(lu_meta == -999.0_dp))
-
-         ! Inside the loop:
+         ! Compute distance using only valid dimensions
          distance = 0.0_dp
          if (use_teff_dim) distance = distance + teff_dist**2
          if (use_logg_dim) distance = distance + logg_dist**2
          if (use_meta_dim) distance = distance + meta_dist**2
 
-         distance = teff_dist**2 + logg_dist**2 + meta_dist**2
-
          do j = 1, 4
             if (distance < min_distances(j)) then
                ! Shift larger distances down
@@ -262,17 +260,17 @@ subroutine interpolate_array(x_in, y_in, x_out, y_out)
       ! Validate input sizes
       if (size(x_in) < 2 .or. size(y_in) < 2) then
          print *, "Error: x_in or y_in arrays have fewer than 2 points."
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       if (size(x_in) /= size(y_in)) then
          print *, "Error: x_in and y_in arrays have different sizes."
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       if (size(x_out) <= 0) then
          print *, "Error: x_out array is empty."
-         stop
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       do i = 1, size(x_out)

colors/private/linear_interp.f90

@@ -24,6 +24,7 @@
 module linear_interp
    use const_def, only: dp
    use colors_utils, only: dilute_flux
+   use utils_lib, only: mesa_error
    implicit none
 
    private
@@ -68,15 +69,17 @@ subroutine construct_sed_linear(teff, log_g, metallicity, R, d, file_names, &
       INQUIRE (file=bin_filename, EXIST=file_exists)
 
       if (.not. file_exists) then
-         stop 'Missing required binary file for interpolation'
+         print *, 'Missing required binary file for interpolation'
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       ! Load the data from binary file
       call load_binary_data(bin_filename, teff_grid, logg_grid, meta_grid, &
                             wavelengths, precomputed_flux_cube, status)
 
       if (status /= 0) then
-         stop 'Binary data loading error'
+         print *, 'Binary data loading error'
+         call mesa_error(__FILE__, __LINE__)
       end if
 
       n_teff = size(teff_grid)