Commit 6782de6f authored by jonas's avatar jonas
Browse files

cfitsio reads different order than astropy

parent 835e34bf
......@@ -76,7 +76,7 @@ def demod_hrt(data,pmp_temp):
def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, bit_flat = True, norm_f = True, clean_f = False,
sigma = 59, flat_states = 24, prefilter_f = None,flat_c = True, dark_c = True, fs_c = True,
demod = True, norm_stokes = True, out_dir = './', out_demod_file = False, out_demod_filename = None,
ItoQUV = False, ctalk_params = None, rte = False, out_rte_filename = None, pmilos = True):
ItoQUV = False, ctalk_params = None, rte = False, out_rte_filename = None, p_milos = True):
'''
PHI-HRT data reduction pipeline
......@@ -148,7 +148,7 @@ def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, bit_flat =
invert using cmilos, options: 'RTE' for Milne Eddington Inversion, 'CE' for Classical Estimates, 'CE+RTE' for combined
out_rte_filename: str, DEFAULT = ''
if '', takes last 10 characters of input scan filename (assumes its a DID), change if want other name
pmilos: bool, DEFAULT = True
p_milos: bool, DEFAULT = True
if True, will execute the RTE inversion using the parallel version of the CMILOS code on 16 processors
Returns
......@@ -876,7 +876,7 @@ def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, bit_flat =
if rte == 'RTE' or rte == 'CE' or rte == 'CE+RTE':
if pmilos:
if p_milos:
try:
pmilos(data_f, wve_axis_arr, data_shape, cpos_arr, data, rte, field_stop, start_row, start_col, out_rte_filename, out_dir)
......
......@@ -50,7 +50,7 @@ out_names = ['0024160031000_test_rte']#, '0024160032000_noflat', '0024160033000_
phihrt_pipe(sciencedata_fits_filenames, flat_f = flatfield_fits_filename, dark_f = darkfield_fits_filename, scale_data = False,
bit_flat = True, norm_f = True, clean_f = False, sigma = 59, flat_states = 24, norm_stokes = True, prefilter_f = None,
dark_c = True, flat_c = True, fs_c = True, demod = True, ctalk_params = c_talk_params, ItoQUV = False, out_demod_file = True,
out_demod_filename = out_names, out_dir = '/data/slam/home/sinjan/hrt_pipe_results/april_2020/', rte = 'RTE', out_rte_filename='', pmilos = False)
out_demod_filename = out_names, out_dir = '/data/slam/home/sinjan/hrt_pipe_results/april_2020_pmilos/', rte = 'RTE', out_rte_filename=None, p_milos = True)
"""
Input Parameters:
----------
......
......@@ -312,13 +312,15 @@ def pmilos(data_f, wve_axis_arr, data_shape, cpos_arr, data, rte, field_stop, st
#write wavelengths to wavelength.fits file for the settings
wave_input = np.zeros((6,2))
wave_input[:,0] = int(1)
wave_input[:,1] = wave_axis
wave_input = np.zeros((2,6)) #cfitsio reads dimensions in opposite order
wave_input[0,:] = 1
wave_input[1,:] = wave_axis
print(wave_axis)
hdr = fits.Header()
primary_hdu = fits.PrimaryH0DU(wave_input, header = hdr)
primary_hdu = fits.PrimaryHDU(wave_input, header = hdr)
hdul = fits.HDUList([primary_hdu])
hdul.writeto(f'./P-MILOS/run/wavelength_tmp.fits', overwrite=True)
......@@ -342,26 +344,32 @@ def pmilos(data_f, wve_axis_arr, data_shape, cpos_arr, data, rte, field_stop, st
cwd = os.getcwd()
os.chdir("./P-MILOS/run/")
cmd = "mpiexec -np 16 ../pmilos.x pmilos.minit" #PMILOS_LOC+"./milos"
cmd = "mpiexec -np 10 ../pmilos.x pmilos.minit" #PMILOS_LOC+"./milos"
if rte == 'RTE':
cmd = "mpiexec -np 16 ../pmilos.x pmilos.minit"
cmd = "mpiexec -np 10 ../pmilos.x pmilos.minit"
if rte == 'CE':
cmd = "mpiexec -np 16 ../pmilos.x pmilos_ce.minit"
cmd = "mpiexec -np 10 ../pmilos.x pmilos_ce.minit"
if rte == 'CE+RTE':
print("CE+RTE not possible on PMILOS, performing RTE instead")
cmd = "mpiexec -np 16 ../pmilos.x pmilos.minit"
cmd = "mpiexec -np 10 ../pmilos.x pmilos.minit"
rte_on = subprocess.call(cmd,shell=True)
os.chdir(cwd)
with fits.open('./P-MILOS/run/results/inv_input_tmp_mod.fits') as hdu_list:
if rte == 'CE':
out_file = 'inv_input_tmp_mod_ce.fits'
else:
out_file = 'inv_input_tmp_mod.fits'
with fits.open(f'./P-MILOS/run/results/{out_file}') as hdu_list:
result = hdu_list[0].data
del_dummy = subprocess.call("rm ./P-MILOS/run/results/inv_input_tmp_mod.fits",shell=True) #must delete the output file
del_dummy = subprocess.call(f"rm ./P-MILOS/run/results/{out_file}.fits",shell=True) #must delete the output file
#result has dimensions [rows,cols,13]
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment