Skip to content

Commit

Permalink
LSF + Fine Time
Browse files Browse the repository at this point in the history
  • Loading branch information
elizabethg60 committed Jan 15, 2025
1 parent ba560e1 commit fa39c78
Show file tree
Hide file tree
Showing 17 changed files with 9,557 additions and 97 deletions.
Binary file added __pycache__/NEID_LSF.cpython-310.pyc
Binary file not shown.
781 changes: 781 additions & 0 deletions data/NEID_October_15s_time.csv

Large diffs are not rendered by default.

7,281 changes: 7,281 additions & 0 deletions data/NEID_October_fine_time.csv

Large diffs are not rendered by default.

867 changes: 867 additions & 0 deletions data/neidMaster_LSFParameterCoefficients20230824HRSci_v1.txt

Large diffs are not rendered by default.

394 changes: 348 additions & 46 deletions scripts/NEID_GRASS_all_lines.jl

Large diffs are not rendered by default.

164 changes: 164 additions & 0 deletions scripts/NEID_LSF.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,164 @@
from astropy.io import fits
import matplotlib.pyplot as plt
import numpy as np
from scipy.special import erf

import sys
import os

# Add the directory containing the script to sys.path
sys.path.append(os.path.abspath(os.path.dirname(__file__)))

def tophat_gauss_kernel( x, mean, fwhm, box_width ):
"""A function to return a top hat convolved Gaussian kernel provided the two width parameters.
Note that the kernel values are not normalized (needed for convolution).
Parameters
----------
x : array
An array of pixel locations to compute the kernel at.
mean : float
The mean of the kernel at the point in the observed spectrum to convolve at.
fwhm : float
The Gaussian FWHM of the kernel at the point in the observed spectrum to convolve at.
box_width : float
The top hat box width of the kernel at the point in the observed spectrum to convolve at.
Returns
-------
kernel_values : array
The values of the LSF convolution kernel.
"""

# First convert the FWHM to a sigma
sigma = fwhm / ( 2 * np.sqrt( 2 * np.log( 2 ) ) )

# Generate the arguments for the two erf components.
erf_arg_1 = ( ( 2 * mean + box_width ) - 2 * x ) / ( 2 * np.sqrt( 2 ) * sigma )
erf_arg_2 = ( (-2 * mean + box_width ) + 2 * x ) / ( 2 * np.sqrt( 2 ) * sigma )

# Generate the kernel values
kernel_values = erf( erf_arg_1 ) + erf( erf_arg_2 )

return kernel_values

def read_lsf_masterfile(file_name):
"""Reads in a text file (in form of IRAF beam trace file) as a dictionary.
Taken from the PRVextract.libs.beamtrace library file in the NEID DRP.
Parameters
----------
file_name : str
The path to the text file to read in.
Returns
-------
file_dict : dict of dicts
The dictionary corresponding to the input file.
"""

# Initialize output dictionary
file_dict = {}

# Go through each of the lines in the file to append to output dictionary
current_aperture = None
for line in open(file_name):
line = line.split()

# Skip empty lines
if not line:
continue

# Read the line
key = line[0]
value = [parseValue(value) for value in line[1:]]

# If only one value for line, turn from list to single value
if len(value) == 1:
value = value[0]

# Insert the information into the dictionary
if key == 'APERTURE':
if current_aperture is not None:
file_dict[current_aperture['APERTURE']] = current_aperture
current_aperture = dict({'APERTURE': value})
else:
current_aperture[key] = value

if current_aperture is not None:
file_dict[current_aperture['APERTURE']] = current_aperture

return file_dict

def parseValue(value):
"""Utility function for parsing a value when reading in beam trace file.
"""

for parse_function in [int, float, str]:
try:
return parse_function(value)
except:
continue

return value

def get_variable_lsf_kernel_values(obs_pixel, kernel_pixel_arr, wavelength_spacing, order):
"""Function to generate a 1D LSF kernel array given parameterization defined in the input masterfile.
The LSF profile has parameters that are fit with polynomials in pixel space across each order.
Parameters
----------
obs_pixel : float
The pixel location within spectral order to generate LSF at.
kernel_pixel_arr : array.
The pixel locations at which the kernel will be evaluated.
kernel_profile_information : dict
Dictionary with GLOBAL information about the LSF profile parameterization.
kernel_parameters : dict
The order-specific dictionary from the input LSF masterfile containing the kernel parameter fit information.
wavelength_spacing : float
The equal spacing of the wavelength grid being convolved, for normalizing the kernel.
Returns
-------
lsf_kernel_values : array
The values of the LSF convolution kernel for the input kernel_pixel_arr.
"""
lsf_info = read_lsf_masterfile("data/neidMaster_LSFParameterCoefficients20230824HRSci_v1.txt")
kernel_profile_information = lsf_info['GLOBAL']
kernel_parameters = lsf_info[order]
# Scale the input observation pixel value to be in the domain over which the LSF parameters are fit (turning [0,9215] to [-1,1])
obs_pixel_m1top1 = scale_interval_m1top1(obs_pixel, *kernel_profile_information['DOMAIN'])

### Pick out the right line profile parameterization from the input master file. Currently just top hat gaussian
if kernel_profile_information['PROFILE'] == 'TOP_HAT_GAUSSIAN':

# Take the polynomial coefficients describing the kernel parameters vs. pixel for this order
kernel_fwhm_poly = kernel_parameters['COEFF_1'] # Coefficient 1 is the Gaussian FWHM
kernel_boxwidth_poly = kernel_parameters['COEFF_2'] # Coefficient 2 is the top hat box width

# Now apply the fits to the input observed pixel value
kernel_fwhm = np.polyval( kernel_fwhm_poly, obs_pixel_m1top1 )
kernel_boxwidth = np.polyval( kernel_boxwidth_poly, obs_pixel_m1top1 )

# Generate the actual kernel values
lsf_kernel_values = tophat_gauss_kernel(kernel_pixel_arr, obs_pixel, kernel_fwhm, kernel_boxwidth)

# Normalize the kernel values so that their integral is 1 (for flux conservation)
lsf_kernel_values /= np.trapz(lsf_kernel_values, dx = wavelength_spacing)

return lsf_kernel_values

def scale_interval_m1top1(x,a,b):
"""Scales input x values over interval a to b onto the range of -1 to 1.
"""

return (2.0 * x - (b + a)) / (b - a)
1 change: 1 addition & 0 deletions scripts/NEID_RV_ccf.jl
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ vacwav = λ_air_to_vac.(airwav)

# line information for identification
orders = [57, 57, 60, 60, 60, 60, 61, 61, 61, 61, 61, 61, 64, 64, 74, 74, 75, 75, 75, 75, 77, 77]
pixels = [3021:3086, 3062:3128, 2227:2287, 2355:2416, 2464:2526, 2564:2626, 2702:2764, 2738:2801, 2881:2944, 3003:3067, 3057:3093, 3070:3134, 2303:2362, 2529:2589, 3563:3622, 3765:3824, 378:425, 414:461, 478:526, 674:722, 734:781, 802:849]

datadir = string(abspath("data"))
df = CSV.read(joinpath(datadir, "optimized_depth.csv"), DataFrame)
Expand Down
88 changes: 57 additions & 31 deletions scripts/debug_gpu_eclipse.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ using SPICE
using CSV
using DataFrames
import PyPlot
using Statistics
plt = PyPlot

#NEID location
Expand All @@ -12,7 +13,7 @@ obs_long = -111.5967
alt = 2.097938

#convert from utc to et as needed by SPICE
time_stamps = utc2et.(["2023-10-14T14:29:23.500000"])
time_stamps = utc2et.(["2023-10-14T15:26:45.500000", "2023-10-14T15:28:07.500000", "2023-10-14T15:29:30.500000", "2023-10-14T15:30:53.500000", "2023-10-14T15:32:15.500000", "2023-10-14T15:33:38.500000", "2023-10-14T15:35:01.500000", "2023-10-14T15:36:23.500000", "2023-10-14T15:37:46.500000", "2023-10-14T15:39:09.500000", "2023-10-14T15:40:31.500000", "2023-10-14T15:41:54.500000", "2023-10-14T15:43:17.500000", "2023-10-14T15:44:39.500000", "2023-10-14T15:46:02.500000", "2023-10-14T15:47:25.500000", "2023-10-14T15:48:47.500000", "2023-10-14T15:50:10.500000", "2023-10-14T15:51:33.500000", "2023-10-14T15:52:56.500000", "2023-10-14T15:54:18.500000", "2023-10-14T15:55:41.500000", "2023-10-14T15:57:04.500000", "2023-10-14T15:58:26.500000", "2023-10-14T15:59:49.500000", "2023-10-14T16:01:12.500000", "2023-10-14T16:02:34.500000", "2023-10-14T16:03:57.500000", "2023-10-14T16:05:20.500000", "2023-10-14T16:06:42.500000", "2023-10-14T16:08:05.500000", "2023-10-14T16:09:28.500000", "2023-10-14T16:10:50.500000", "2023-10-14T16:12:13.500000", "2023-10-14T16:13:36.500000", "2023-10-14T16:14:58.500000", "2023-10-14T16:16:21.500000", "2023-10-14T16:17:44.500000", "2023-10-14T16:19:06.500000", "2023-10-14T16:20:29.500000", "2023-10-14T16:21:52.500000", "2023-10-14T16:23:15.500000", "2023-10-14T16:24:37.500000", "2023-10-14T16:26:00.500000", "2023-10-14T16:27:23.500000", "2023-10-14T16:28:45.500000", "2023-10-14T16:30:08.500000", "2023-10-14T16:31:31.500000", "2023-10-14T16:32:53.500000", "2023-10-14T16:34:16.500000", "2023-10-14T16:35:39.500000", "2023-10-14T16:37:01.500000", "2023-10-14T16:38:24.500000", "2023-10-14T16:39:47.500000", "2023-10-14T16:41:09.500000", "2023-10-14T16:42:32.500000", "2023-10-14T16:43:55.500000", "2023-10-14T16:45:17.500000", "2023-10-14T16:46:40.500000", "2023-10-14T16:48:03.500000", "2023-10-14T16:49:25.500000", "2023-10-14T16:50:48.500000", "2023-10-14T16:52:11.500000", "2023-10-14T16:53:33.500000", "2023-10-14T16:54:56.500000", "2023-10-14T16:56:19.500000", "2023-10-14T16:57:42.500000", "2023-10-14T16:59:04.500000", "2023-10-14T17:00:27.500000", "2023-10-14T17:01:50.500000", "2023-10-14T17:03:12.500000", "2023-10-14T17:04:35.500000", "2023-10-14T17:05:58.500000", "2023-10-14T17:07:20.500000", "2023-10-14T17:08:43.500000", "2023-10-14T17:10:06.500000", "2023-10-14T17:11:28.500000", "2023-10-14T17:12:51.500000", "2023-10-14T17:14:14.500000", "2023-10-14T17:15:36.500000", "2023-10-14T17:16:59.500000", "2023-10-14T17:18:22.500000", "2023-10-14T17:19:44.500000", "2023-10-14T17:21:07.500000", "2023-10-14T17:22:30.500000", "2023-10-14T17:23:52.500000", "2023-10-14T17:25:15.500000", "2023-10-14T17:26:38.500000", "2023-10-14T17:28:01.500000", "2023-10-14T17:29:23.500000", "2023-10-14T17:30:46.500000", "2023-10-14T17:32:09.500000", "2023-10-14T17:33:31.500000", "2023-10-14T17:34:54.500000", "2023-10-14T17:36:17.500000", "2023-10-14T17:37:39.500000", "2023-10-14T17:39:02.500000", "2023-10-14T17:40:25.500000", "2023-10-14T17:41:47.500000", "2023-10-14T17:43:10.500000", "2023-10-14T17:44:33.500000", "2023-10-14T17:45:55.500000", "2023-10-14T17:47:18.500000", "2023-10-14T17:48:41.500000", "2023-10-14T17:50:03.500000", "2023-10-14T17:51:26.500000", "2023-10-14T17:52:49.500000", "2023-10-14T17:54:11.500000", "2023-10-14T17:55:34.500000", "2023-10-14T17:56:57.500000", "2023-10-14T17:58:20.500000", "2023-10-14T17:59:42.500000", "2023-10-14T18:01:05.500000", "2023-10-14T18:02:28.500000", "2023-10-14T18:03:50.500000", "2023-10-14T18:05:13.500000", "2023-10-14T18:06:36.500000", "2023-10-14T18:07:58.500000", "2023-10-14T18:09:21.500000", "2023-10-14T18:10:44.500000", "2023-10-14T18:12:06.500000", "2023-10-14T18:13:29.500000", "2023-10-14T18:14:52.500000", "2023-10-14T18:16:14.500000", "2023-10-14T18:17:37.500000", "2023-10-14T18:19:00.500000", "2023-10-14T18:20:22.500000", "2023-10-14T18:21:45.500000", "2023-10-14T18:23:08.500000", "2023-10-14T18:24:30.500000", "2023-10-14T18:25:53.500000", "2023-10-14T18:27:16.500000", "2023-10-14T18:28:38.500000", "2023-10-14T18:30:01.500000", "2023-10-14T18:31:24.500000", "2023-10-14T18:32:47.500000", "2023-10-14T18:34:09.500000", "2023-10-14T18:35:32.500000", "2023-10-14T18:36:55.500000", "2023-10-14T18:38:17.500000", "2023-10-14T18:39:40.500000", "2023-10-14T18:41:03.500000", "2023-10-14T18:42:25.500000", "2023-10-14T18:43:48.500000", "2023-10-14T18:45:11.500000", "2023-10-14T18:46:33.500000", "2023-10-14T18:47:56.500000", "2023-10-14T18:49:19.500000", "2023-10-14T18:50:41.500000", "2023-10-14T18:52:04.500000", "2023-10-14T18:53:27.500000", "2023-10-14T18:54:49.500000", "2023-10-14T18:56:12.500000", "2023-10-14T18:57:35.500000", "2023-10-14T18:58:57.500000", "2023-10-14T19:00:20.500000", "2023-10-14T19:01:43.500000", "2023-10-14T19:03:06.500000"])

N = 50
Nt = length(time_stamps)
Expand Down Expand Up @@ -46,34 +47,59 @@ z_rot_matrix = Matrix{Matrix{Float64}}(undef, length(disk.ϕc), maximum(disk.Nθ
zenith_matrix = Matrix{Matrix{Float64}}(undef, length(disk.ϕc), maximum(disk.Nθ))

# compute geometry for timestamp
t = 1
for t in 1:length(time_stamps)

GRASS.eclipse_compute_quantities!(disk, time_stamps[t], t, obs_long, obs_lat, alt, wsp.ϕc, wsp.θc,
wsp.μs, wsp.z_rot, wsp.dA, wsp.xyz, wsp.ax_codes,
mem.dA_total_proj_mean, mem.mean_weight_v_no_cb,
mem.mean_weight_v_earth_orb, mem.pole_vector_grid,
mem.SP_sun_pos, mem.SP_sun_vel, mem.SP_bary, mem.SP_bary_pos,
mem.SP_bary_vel, mem.OP_bary, mem.mu_grid, mem.projected_velocities_no_cb,
mem.distance, mem.v_scalar_grid, mem.v_earth_orb_proj, mu_grid_matrix,
dA_total_proj_matrix, idx1_matrix, idx3_matrix, z_rot_matrix, zenith_matrix)

lambdas_cpu, outspec_cpu = GRASS.eclipse_compute_intensity(disk, lines, 1.0, LD_type, idx1_matrix, idx3_matrix,
mu_grid_matrix, z_rot_matrix, dA_total_proj_matrix, wsp.ld, wsp.z_rot, zenith_matrix,
wsp.μs, wsp.ax_codes, wsp.dA, wsp.μs, wsp.ax_codes, wsp.dA, true, t, wsp.ext)

"""
GPU GEOMETRY CODE CALLS BELOW
"""
gpu_allocs = GRASS.GPUAllocsEclipse(spec, disk, 1)

GRASS.calc_eclipse_quantities_gpu!(time_stamps[t], obs_long, obs_lat, alt, lines, LD_type, 1.0, 1.0, disk, gpu_allocs)

# make a plot
plt.imshow(wsp.ext .- Array(gpu_allocs.ext))
# plt.imshow(wsp.z_rot .* GRASS.c_ms, vmin=-2000, vmax=2000, cmap="seismic")
# plt.imshow(Array(gpu_allocs.z_rot) .* GRASS.c_ms .- wsp.z_rot .* GRASS.c_ms)
plt.title("ext")
plt.colorbar()
plt.savefig("gpu_test.png")
plt.show()
plt.clf()
GRASS.eclipse_compute_quantities!(disk, time_stamps[t], t, obs_long, obs_lat, alt, wsp.ϕc, wsp.θc,
wsp.μs, wsp.z_rot, wsp.dA, wsp.xyz, wsp.ax_codes,
mem.dA_total_proj_mean, mem.mean_weight_v_no_cb,
mem.mean_weight_v_earth_orb, mem.pole_vector_grid,
mem.SP_sun_pos, mem.SP_sun_vel, mem.SP_bary, mem.SP_bary_pos,
mem.SP_bary_vel, mem.OP_bary, mem.mu_grid, mem.projected_velocities_no_cb,
mem.distance, mem.v_scalar_grid, mem.v_earth_orb_proj, mu_grid_matrix,
dA_total_proj_matrix, idx1_matrix, idx3_matrix, z_rot_matrix, zenith_matrix)

# lambdas_cpu, outspec_cpu = GRASS.eclipse_compute_intensity(disk, lines, 1.0, LD_type, idx1_matrix, idx3_matrix,
# mu_grid_matrix, z_rot_matrix, dA_total_proj_matrix, wsp.ld, wsp.z_rot, zenith_matrix,
# wsp.μs, wsp.ax_codes, wsp.dA, wsp.μs, wsp.ax_codes, wsp.dA, true, t, wsp.ext)

mean_intensity = GRASS.eclipse_compute_intensity(disk, lines, 1.0, LD_type, idx1_matrix, idx3_matrix,
mu_grid_matrix, mem.mean_weight_v_no_cb[:, :, t], mem.mean_weight_v_earth_orb[:, :, t],
z_rot_matrix, dA_total_proj_matrix, wsp.ld, wsp.z_rot, zenith_matrix, true, wsp.ext)
idx_grid = mean_intensity[:, :, 1] .> 0.0
brightness = mean_intensity[:, :, 1] .* mem.dA_total_proj_mean[:, :, t] .* wsp.ext[:, :, 1]
cheapflux = sum(view(brightness, idx_grid))

final_weight_v_no_cb = sum(view(mem.mean_weight_v_no_cb[:, :, t] .* brightness, idx_grid)) / cheapflux
final_weight_v_no_cb += mean(view(mem.mean_weight_v_earth_orb[:, :, t], idx_grid))

"""
GPU GEOMETRY CODE CALLS BELOW
"""
gpu_allocs = GRASS.GPUAllocsEclipse(spec, disk, 1)

GRASS.calc_eclipse_quantities_gpu!(time_stamps[t], obs_long, obs_lat, alt, lines, LD_type, 1.0, 1.0, disk, gpu_allocs)

idx_grid2 = Array(gpu_allocs.ld[:, :, 1]) .> 0.0

brightness2 = Array(gpu_allocs.ld[:, :, 1]) .* Array(gpu_allocs.dA) .* Array(gpu_allocs.ext[:, :, 1])

cheapflux2 = sum(view(brightness2, idx_grid2))

final_weight_v_no_cb2 = sum(view(Array(gpu_allocs.projected_v) .* brightness2, idx_grid2)) / cheapflux2
final_weight_v_no_cb2 += mean(view(Array(gpu_allocs.earth_v), idx_grid2))

# plt.scatter(t, final_weight_v_no_cb - final_weight_v_no_cb2, color = "r")
println(final_weight_v_no_cb - final_weight_v_no_cb2)
end
# plt.savefig("test")
# plt.clf()

# # make a plot
# plt.imshow(mem.mean_weight_v_earth_orb[:, :, t] .- Array(gpu_allocs.earth_v))
# # plt.imshow(wsp.z_rot .* GRASS.c_ms, vmin=-2000, vmax=2000, cmap="seismic")
# # plt.imshow(Array(gpu_allocs.z_rot) .* GRASS.c_ms .- wsp.z_rot .* GRASS.c_ms)
# plt.title("earth_v")
# plt.colorbar()
# plt.savefig("gpu_test.png")
# plt.show()
# plt.clf()
2 changes: 1 addition & 1 deletion projected_RV.jl → scripts/projected_RV.jl
Original file line number Diff line number Diff line change
Expand Up @@ -153,12 +153,12 @@ function projected_RV_gpu(time, LD_type, ext_toggle)
# make the spec composite type instances
spec = GRASS.SpecParams(lines=lines, depths=depths, variability=variability,
blueshifts=blueshifts, templates=templates, resolution=resolution)
gpu_allocs = GRASS.GPUAllocsEclipse(spec, disk, Int(length(lines)), precision=Float64, verbose=true)

RV_list_no_cb = Vector{Float64}(undef,size(time_stamps)...)
intensity_list = Vector{Float64}(undef,size(time_stamps)...)
mean_intensity_list = Vector{Matrix{Float64}}(undef,size(time_stamps)...)
for t in 1:disk.Nt
gpu_allocs = GRASS.GPUAllocsEclipse(spec, disk, Int(length(lines)), precision=Float64, verbose=true)
if t < 25
neid_ext_coeff = extinction_coeff[extinction_coeff[!, "Wavelength"] .== λrest[i], "Ext1"][1]
elseif t >= 25 && t < 46
Expand Down
4 changes: 2 additions & 2 deletions scripts/slurm_gpu.slurm
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,9 @@
#SBATCH --gpus-per-task=1
#SBATCH --mem-per-cpu=64GB
#SBATCH --time=24:00:00
#SBATCH --job-name=N_convergence
#SBATCH --job-name=FineTime
#SBATCH --chdir=/storage/home/efg5335/work/GRASS
#SBATCH --output=/storage/home/efg5335/work/GRASS/N_convergence.out
#SBATCH --output=/storage/home/efg5335/work/GRASS/FineTime.out

echo "About to start: $SLURM_JOB_NAME"
date
Expand Down
5 changes: 3 additions & 2 deletions src/convenience_eclipse.jl
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,7 @@ function synth_Eclipse_gpu(spec::SpecParams{T}, disk::DiskParamsEclipse{T},
# allocate memory needed for rossiter computations
gpu_allocs = GPUAllocsEclipse(spec, disk, Int(length(wavelength)), precision=precision, verbose=verbose)

brightness_arr = Vector{Vector{Float64}}(undef,size(templates)...)
# run the simulation and return
for (idx, file) in enumerate(templates)
# get temporary specparams with lines for this run
Expand All @@ -117,9 +118,9 @@ function synth_Eclipse_gpu(spec::SpecParams{T}, disk::DiskParamsEclipse{T},
soldata = GPUSolarData(soldata_cpu, precision=precision)

# run the simulation and multiply flux by this spectrum
disk_sim_eclipse_gpu(spec_temp, disk, soldata, gpu_allocs, flux,
brightness_arr[idx] = disk_sim_eclipse_gpu(spec_temp, disk, soldata, gpu_allocs, flux,
obs_long, obs_lat, alt, time_stamps, wavelength,
ext_coeff, ext_toggle, LD_type, skip_times=skip_times)
end
return spec.lambdas, flux
return spec.lambdas, flux, brightness_arr
end
4 changes: 2 additions & 2 deletions src/disk_sim_eclipse.jl
Original file line number Diff line number Diff line change
Expand Up @@ -101,13 +101,13 @@ function disk_sim_eclipse(spec::SpecParams{T}, disk::DiskParamsEclipse{T}, solda
wsp.widt .= copy(view(soldata.wid[key], :, tloop[i,j]))

# get amount of convective blueshift needed
# extra_z = spec.conv_blueshifts[l] - z_cbs_avg
extra_z = spec.conv_blueshifts[l] - z_cbs_avg

# get shifted line center
λΔD = spec.lines[l]
λΔD *= (1.0 + z_rot)
λΔD *= (1.0 + z_cbs)
# λΔD *= (1.0 + extra_z)
λΔD *= (1.0 + extra_z)

# fix bisector and width if variability is turned off
if !spec.variability[l]
Expand Down
Loading

0 comments on commit fa39c78

Please sign in to comment.