Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
jonas.sinjan
hrt_pipeline
Commits
e70e186e
Commit
e70e186e
authored
Jul 21, 2021
by
jonas
Browse files
pmilos init
parent
ab211a89
Changes
5
Hide whitespace changes
Inline
Side-by-side
.gitignore
View file @
e70e186e
...
...
@@ -15,4 +15,5 @@ dummy*.txt
_RTE.npz
*.npz
.env
run.py
\ No newline at end of file
run.py
P-MILOS/*
\ No newline at end of file
README.md
View file @
e70e186e
...
...
@@ -64,7 +64,7 @@ OR : use `download_files.py` to download images from the attic repository: https
1.
Compile milos:
```
bash
make clea
n
make clea
r
make
```
...
...
hrt_pipe.py
View file @
e70e186e
...
...
@@ -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
=
''
):
ItoQUV
=
False
,
ctalk_params
=
None
,
rte
=
False
,
out_rte_filename
=
None
,
pmilos
=
True
):
'''
PHI-HRT data reduction pipeline
...
...
@@ -96,7 +96,7 @@ def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, bit_flat =
13. calibration
a) ghost correction - not implemented yet
b) cross talk correction
14. rte inversion with cmilos (CE, RTE or CE+RTE)
14. rte inversion with
pmilos/
cmilos (CE, RTE or CE+RTE)
a) output rte data products to fits file
Parameters
...
...
@@ -148,6 +148,8 @@ 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
if True, will execute the RTE inversion using the parallel version of the CMILOS code on 16 processors
Returns
-------
...
...
@@ -874,162 +876,17 @@ def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, bit_flat =
if
rte
==
'RTE'
or
rte
==
'CE'
or
rte
==
'CE+RTE'
:
print
(
" "
)
printc
(
'-->>>>>>> RUNNING CMILOS '
,
color
=
bcolors
.
OKGREEN
)
try
:
CMILOS_LOC
=
os
.
path
.
realpath
(
__file__
)
CMILOS_LOC
=
CMILOS_LOC
[:
-
11
]
+
'cmilos/'
#11 as hrt_pipe.py is 11 characters
if
os
.
path
.
isfile
(
CMILOS_LOC
+
'milos'
):
printc
(
"Cmilos executable located at:"
,
CMILOS_LOC
,
color
=
bcolors
.
WARNING
)
else
:
raise
ValueError
(
'Cannot find cmilos:'
,
CMILOS_LOC
)
except
ValueError
as
err
:
printc
(
err
.
args
[
0
],
color
=
bcolors
.
FAIL
)
printc
(
err
.
args
[
1
],
color
=
bcolors
.
FAIL
)
return
wavelength
=
6173.3356
for
scan
in
range
(
int
(
data_shape
[
-
1
])):
start_time
=
time
.
time
()
if
isinstance
(
data_f
,
str
):
file_path
=
data_f
elif
isinstance
(
data_f
,
list
):
file_path
=
data_f
[
scan
]
wave_axis
=
wve_axis_arr
[
scan
]
#must invert each scan independently, as cmilos only takes in one dataset at a time
#get wave_axis from the header information of the science scans
if
cpos_arr
[
0
]
==
0
:
shift_w
=
wave_axis
[
3
]
-
wavelength
elif
cpos_arr
[
0
]
==
5
:
shift_w
=
wave_axis
[
2
]
-
wavelength
wave_axis
=
wave_axis
-
shift_w
print
(
'It is assumed the wavelength array is given by the header'
)
#print(wave_axis,color = bcolors.WARNING)
print
(
"Wave axis is: "
,
(
wave_axis
-
wavelength
)
*
1000.
)
print
(
'Saving data into dummy_in.txt for RTE input'
)
sdata
=
data
[:,:,:,:,
scan
]
y
,
x
,
p
,
l
=
sdata
.
shape
#print(y,x,p,l)
filename
=
'dummy_in.txt'
with
open
(
filename
,
"w"
)
as
f
:
for
i
in
range
(
x
):
for
j
in
range
(
y
):
for
k
in
range
(
l
):
f
.
write
(
'%e %e %e %e %e
\n
'
%
(
wave_axis
[
k
],
sdata
[
j
,
i
,
0
,
k
],
sdata
[
j
,
i
,
1
,
k
],
sdata
[
j
,
i
,
2
,
k
],
sdata
[
j
,
i
,
3
,
k
]))
#wv, I, Q, U, V
del
sdata
printc
(
f
' ---- >>>>> Inverting data scan number:
{
scan
}
.... '
,
color
=
bcolors
.
OKGREEN
)
cmd
=
CMILOS_LOC
+
"./milos"
cmd
=
fix_path
(
cmd
)
if
rte
==
'RTE'
:
rte_on
=
subprocess
.
call
(
cmd
+
" 6 15 0 0 dummy_in.txt > dummy_out.txt"
,
shell
=
True
)
if
rte
==
'CE'
:
rte_on
=
subprocess
.
call
(
cmd
+
" 6 15 2 0 dummy_in.txt > dummy_out.txt"
,
shell
=
True
)
if
rte
==
'CE+RTE'
:
rte_on
=
subprocess
.
call
(
cmd
+
" 6 15 1 0 dummy_in.txt > dummy_out.txt"
,
shell
=
True
)
#print(rte_on)
printc
(
' ---- >>>>> Reading results.... '
,
color
=
bcolors
.
OKGREEN
)
del_dummy
=
subprocess
.
call
(
"rm dummy_in.txt"
,
shell
=
True
)
#print(del_dummy)
res
=
np
.
loadtxt
(
'dummy_out.txt'
)
npixels
=
res
.
shape
[
0
]
/
12.
#print(npixels)
#print(npixels/x)
result
=
np
.
zeros
((
12
,
y
*
x
)).
astype
(
float
)
rte_invs
=
np
.
zeros
((
12
,
y
,
x
)).
astype
(
float
)
for
i
in
range
(
y
*
x
):
result
[:,
i
]
=
res
[
i
*
12
:(
i
+
1
)
*
12
]
result
=
result
.
reshape
(
12
,
y
,
x
)
result
=
np
.
einsum
(
'ijk->ikj'
,
result
)
rte_invs
=
result
del
result
rte_invs_noth
=
np
.
copy
(
rte_invs
)
"""
From 0 to 11
Counter (PX Id)
Iterations
Strength
Inclination
Azimuth
Eta0 parameter
Doppler width
Damping
Los velocity
Constant source function
Slope source function
Minimum chisqr value
"""
noise_in_V
=
np
.
mean
(
data
[:,:,
3
,
cpos_arr
[
0
],:])
low_values_flags
=
np
.
max
(
np
.
abs
(
data
[:,:,
3
,:,
scan
]),
axis
=-
1
)
<
noise_in_V
# Where values are low
rte_invs
[
2
,
low_values_flags
]
=
0
rte_invs
[
3
,
low_values_flags
]
=
0
rte_invs
[
4
,
low_values_flags
]
=
0
#np.savez_compressed(out_dir+'_RTE', rte_invs=rte_invs, rte_invs_noth=rte_invs_noth)
del_dummy
=
subprocess
.
call
(
"rm dummy_out.txt"
,
shell
=
True
)
#print(del_dummy)
rte_data_products
=
np
.
zeros
((
6
,
rte_invs_noth
.
shape
[
1
],
rte_invs_noth
.
shape
[
1
]))
if
pmilos
:
rte_data_products
[
0
,:,:]
=
rte_invs_noth
[
9
,:,:]
+
rte_invs_noth
[
10
,:,:]
#continuum
rte_data_products
[
1
,:,:]
=
rte_invs_noth
[
2
,:,:]
#b mag strength
rte_data_products
[
2
,:,:]
=
rte_invs_noth
[
3
,:,:]
#inclination
rte_data_products
[
3
,:,:]
=
rte_invs_noth
[
4
,:,:]
#azimuth
rte_data_products
[
4
,:,:]
=
rte_invs_noth
[
8
,:,:]
#vlos
rte_data_products
[
5
,:,:]
=
rte_invs_noth
[
2
,:,:]
*
np
.
cos
(
rte_invs_noth
[
3
,:,:]
*
np
.
pi
/
180.
)
#blos
rte_data_products
*=
field_stop
[
np
.
newaxis
,
start_row
:
start_row
+
data_size
[
0
],
start_col
:
start_col
+
data_size
[
1
]]
#field stop, set outside to 0
if
out_rte_filename
==
''
:
filename_root
=
str
(
file_path
.
split
(
'.fits'
)[
0
][
-
10
:])
else
:
filename_root
=
out_rte_filename
with
fits
.
open
(
file_path
)
as
hdu_list
:
hdu_list
[
0
].
data
=
rte_data_products
hdu_list
.
writeto
(
out_dir
+
filename_root
+
'_rte_data_products.fits'
,
overwrite
=
True
)
with
fits
.
open
(
file_path
)
as
hdu_list
:
hdu_list
[
0
].
data
=
rte_data_products
[
5
,:,:]
hdu_list
.
writeto
(
out_dir
+
filename_root
+
'_blos_rte.fits'
,
overwrite
=
True
)
with
fits
.
open
(
file_path
)
as
hdu_list
:
hdu_list
[
0
].
data
=
rte_data_products
[
4
,:,:]
hdu_list
.
writeto
(
out_dir
+
filename_root
+
'_vlos_rte.fits'
,
overwrite
=
True
)
with
fits
.
open
(
file_path
)
as
hdu_list
:
hdu_list
[
0
].
data
=
rte_data_products
[
0
,:,:]
hdu_list
.
writeto
(
out_dir
+
filename_root
+
'_Icont_rte.fits'
,
overwrite
=
True
)
try
:
pmilos
(
data_f
,
wve_axis_arr
,
data_shape
,
cpos_arr
,
data
,
rte
,
field_stop
,
start_row
,
start_col
,
out_rte_filename
,
out_dir
)
except
ValueError
:
print
(
"Running CMILOS instead!"
)
cmilos
(
data_f
,
wve_axis_arr
,
data_shape
,
cpos_arr
,
data
,
rte
,
field_stop
,
start_row
,
start_col
,
out_rte_filename
,
out_dir
)
printc
(
'--------------------------------------------------------------'
,
bcolors
.
OKGREEN
)
printc
(
f
"------------- CMILOS RTE Run Time:
{
np
.
round
(
time
.
time
()
-
start_time
,
3
)
}
seconds "
,
bcolors
.
OKGREEN
)
printc
(
'--------------------------------------------------------------'
,
bcolors
.
OKGREEN
)
else
:
cmilos
(
data_f
,
wve_axis_arr
,
data_shape
,
cpos_arr
,
data
,
rte
,
field_stop
,
start_row
,
start_col
,
out_rte_filename
,
out_dir
)
else
:
print
(
" "
)
...
...
run.py
View file @
e70e186e
from
hrt_pipe
import
phihrt_pipe
import
numpy
as
np
sciencedata_fits_filenames
=
[
'solo_L0_phi-hrt-ilam_20200420T142022_V202004221451C_0024160031000.fits'
,
'solo_L0_phi-hrt-ilam_20200420T142252_V202004221452C_0024160032000.fits'
,
'solo_L0_phi-hrt-ilam_20200420T142522_V202004221457C_0024160033000.fits'
,
'solo_L0_phi-hrt-ilam_20200420T142752_V202004221511C_0024160034000.fits'
,
'solo_L0_phi-hrt-ilam_20200420T143023_V202004221517C_0024160035000.fits'
,
'solo_L0_phi-hrt-ilam_20200420T143253_V202004221518C_0024160036000.fits'
,
'solo_L0_phi-hrt-ilam_20200420T143523_V202004221522C_0024160037000.fits'
,
'solo_L0_phi-hrt-ilam_20200420T143753_V202004231605C_0024160038000.fits'
,
'solo_L0_phi-hrt-ilam_20200420T144023_V202004231605C_0024160039000.fits'
]
#'/data/slam/home/sinjan/fits_files/solo_L0_phi-hrt-ilam_20200420T141752_V202004221450C_0024160030000.fits'#solo_L1_phi-hrt-ilam_20201117T170209_V202107090920C_0051170001.fits'#solo_L1_phi-hrt-ilam_20201117T170209_V202107060746C_0051170001.fits'#['solo_L1_phi-hrt-ilam_20200528T171109_V202106111600C_0045140102.fits']
sciencedata_fits_filenames
=
[
'solo_L0_phi-hrt-ilam_20200420T142022_V202004221451C_0024160031000.fits'
]
#
,'solo_L0_phi-hrt-ilam_20200420T142252_V202004221452C_0024160032000.fits','solo_L0_phi-hrt-ilam_20200420T142522_V202004221457C_0024160033000.fits','solo_L0_phi-hrt-ilam_20200420T142752_V202004221511C_0024160034000.fits','solo_L0_phi-hrt-ilam_20200420T143023_V202004221517C_0024160035000.fits', 'solo_L0_phi-hrt-ilam_20200420T143253_V202004221518C_0024160036000.fits','solo_L0_phi-hrt-ilam_20200420T143523_V202004221522C_0024160037000.fits', 'solo_L0_phi-hrt-ilam_20200420T143753_V202004231605C_0024160038000.fits','solo_L0_phi-hrt-ilam_20200420T144023_V202004231605C_0024160039000.fits']#'/data/slam/home/sinjan/fits_files/solo_L0_phi-hrt-ilam_20200420T141752_V202004221450C_0024160030000.fits'#solo_L1_phi-hrt-ilam_20201117T170209_V202107090920C_0051170001.fits'#solo_L1_phi-hrt-ilam_20201117T170209_V202107060746C_0051170001.fits'#['solo_L1_phi-hrt-ilam_20200528T171109_V202106111600C_0045140102.fits']
# ['solo_L0_phi-hrt-ilam_20210421T120003_V202106080929C_0144210101.fits', 'solo_L0_phi-hrt-ilam_20210424T120003_V202106141014C_0144240101.fits',
# 'solo_L0_phi-hrt-ilam_20210425T120002_V202106141020C_0144250101.fits', 'solo_L0_phi-hrt-ilam_20210426T120002_V202106162118C_0144260101.fits',
...
...
@@ -45,12 +45,12 @@ c_talk_params[1,0] = q_int
c_talk_params
[
1
,
1
]
=
u_int
c_talk_params
[
1
,
2
]
=
v_int
out_names
=
[
'0024160031000_
noflat'
,
'0024160032000_noflat'
,
'0024160033000_noflat'
,
'0024160034000_noflat'
,
'0024160035000_noflat'
,
'0024160036000_noflat'
,
'0024160037000_noflat'
,
'0024160038000_noflat'
,
'0024160039000_noflat'
]
out_names
=
[
'0024160031000_
test_rte'
]
#
, '0024160032000_noflat', '0024160033000_noflat', '0024160034000_noflat', '0024160035000_noflat', '0024160036000_noflat', '0024160037000_noflat', '0024160038000_noflat', '0024160039000_noflat']
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
=
Fals
e
,
flat_c
=
Fals
e
,
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
=
'
False
'
,
out_rte_filename
=
''
)
dark_c
=
Tru
e
,
flat_c
=
Tru
e
,
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
)
"""
Input Parameters:
----------
...
...
utils.py
View file @
e70e186e
from
astropy.io
import
fits
import
numpy
as
np
import
os
import
time
import
subprocess
class
bcolors
:
HEADER
=
'
\033
[95m'
...
...
@@ -109,4 +112,311 @@ def fix_path(path,dir='forward',verbose=False):
print
(
path
)
return
path
else
:
pass
\ No newline at end of file
pass
def
cmilos
(
data_f
,
wve_axis_arr
,
data_shape
,
cpos_arr
,
data
,
rte
,
field_stop
,
start_row
,
start_col
,
out_rte_filename
,
out_dir
):
print
(
" "
)
printc
(
'-->>>>>>> RUNNING CMILOS '
,
color
=
bcolors
.
OKGREEN
)
try
:
CMILOS_LOC
=
os
.
path
.
realpath
(
__file__
)
CMILOS_LOC
=
CMILOS_LOC
[:
-
8
]
+
'cmilos/'
#-11 as hrt_pipe.py is 11 characters
if
os
.
path
.
isfile
(
CMILOS_LOC
+
'milos'
):
printc
(
"Cmilos executable located at:"
,
CMILOS_LOC
,
color
=
bcolors
.
WARNING
)
else
:
raise
ValueError
(
'Cannot find cmilos:'
,
CMILOS_LOC
)
except
ValueError
as
err
:
printc
(
err
.
args
[
0
],
color
=
bcolors
.
FAIL
)
printc
(
err
.
args
[
1
],
color
=
bcolors
.
FAIL
)
return
wavelength
=
6173.3356
for
scan
in
range
(
int
(
data_shape
[
-
1
])):
start_time
=
time
.
time
()
file_path
=
data_f
[
scan
]
wave_axis
=
wve_axis_arr
[
scan
]
#must invert each scan independently, as cmilos only takes in one dataset at a time
#get wave_axis from the header information of the science scans
if
cpos_arr
[
0
]
==
0
:
shift_w
=
wave_axis
[
3
]
-
wavelength
elif
cpos_arr
[
0
]
==
5
:
shift_w
=
wave_axis
[
2
]
-
wavelength
wave_axis
=
wave_axis
-
shift_w
print
(
'It is assumed the wavelength array is given by the header'
)
#print(wave_axis,color = bcolors.WARNING)
print
(
"Wave axis is: "
,
(
wave_axis
-
wavelength
)
*
1000.
)
print
(
'Saving data into dummy_in.txt for RTE input'
)
sdata
=
data
[:,:,:,:,
scan
]
y
,
x
,
p
,
l
=
sdata
.
shape
#print(y,x,p,l)
filename
=
'dummy_in.txt'
with
open
(
filename
,
"w"
)
as
f
:
for
i
in
range
(
x
):
for
j
in
range
(
y
):
for
k
in
range
(
l
):
f
.
write
(
'%e %e %e %e %e
\n
'
%
(
wave_axis
[
k
],
sdata
[
j
,
i
,
0
,
k
],
sdata
[
j
,
i
,
1
,
k
],
sdata
[
j
,
i
,
2
,
k
],
sdata
[
j
,
i
,
3
,
k
]))
#wv, I, Q, U, V
del
sdata
printc
(
f
' ---- >>>>> Inverting data scan number:
{
scan
}
.... '
,
color
=
bcolors
.
OKGREEN
)
cmd
=
CMILOS_LOC
+
"./milos"
cmd
=
fix_path
(
cmd
)
if
rte
==
'RTE'
:
rte_on
=
subprocess
.
call
(
cmd
+
" 6 15 0 0 dummy_in.txt > dummy_out.txt"
,
shell
=
True
)
if
rte
==
'CE'
:
rte_on
=
subprocess
.
call
(
cmd
+
" 6 15 2 0 dummy_in.txt > dummy_out.txt"
,
shell
=
True
)
if
rte
==
'CE+RTE'
:
rte_on
=
subprocess
.
call
(
cmd
+
" 6 15 1 0 dummy_in.txt > dummy_out.txt"
,
shell
=
True
)
#print(rte_on)
printc
(
' ---- >>>>> Reading results.... '
,
color
=
bcolors
.
OKGREEN
)
del_dummy
=
subprocess
.
call
(
"rm dummy_in.txt"
,
shell
=
True
)
#print(del_dummy)
res
=
np
.
loadtxt
(
'dummy_out.txt'
)
npixels
=
res
.
shape
[
0
]
/
12.
#print(npixels)
#print(npixels/x)
result
=
np
.
zeros
((
12
,
y
*
x
)).
astype
(
float
)
rte_invs
=
np
.
zeros
((
12
,
y
,
x
)).
astype
(
float
)
for
i
in
range
(
y
*
x
):
result
[:,
i
]
=
res
[
i
*
12
:(
i
+
1
)
*
12
]
result
=
result
.
reshape
(
12
,
y
,
x
)
result
=
np
.
einsum
(
'ijk->ikj'
,
result
)
rte_invs
=
result
del
result
rte_invs_noth
=
np
.
copy
(
rte_invs
)
"""
From 0 to 11
Counter (PX Id)
Iterations
Strength
Inclination
Azimuth
Eta0 parameter
Doppler width
Damping
Los velocity
Constant source function
Slope source function
Minimum chisqr value
"""
noise_in_V
=
np
.
mean
(
data
[:,:,
3
,
cpos_arr
[
0
],:])
low_values_flags
=
np
.
max
(
np
.
abs
(
data
[:,:,
3
,:,
scan
]),
axis
=-
1
)
<
noise_in_V
# Where values are low
rte_invs
[
2
,
low_values_flags
]
=
0
rte_invs
[
3
,
low_values_flags
]
=
0
rte_invs
[
4
,
low_values_flags
]
=
0
#np.savez_compressed(out_dir+'_RTE', rte_invs=rte_invs, rte_invs_noth=rte_invs_noth)
del_dummy
=
subprocess
.
call
(
"rm dummy_out.txt"
,
shell
=
True
)
#print(del_dummy)
rte_data_products
=
np
.
zeros
((
6
,
rte_invs_noth
.
shape
[
1
],
rte_invs_noth
.
shape
[
2
]))
rte_data_products
[
0
,:,:]
=
rte_invs_noth
[
9
,:,:]
+
rte_invs_noth
[
10
,:,:]
#continuum
rte_data_products
[
1
,:,:]
=
rte_invs_noth
[
2
,:,:]
#b mag strength
rte_data_products
[
2
,:,:]
=
rte_invs_noth
[
3
,:,:]
#inclination
rte_data_products
[
3
,:,:]
=
rte_invs_noth
[
4
,:,:]
#azimuth
rte_data_products
[
4
,:,:]
=
rte_invs_noth
[
8
,:,:]
#vlos
rte_data_products
[
5
,:,:]
=
rte_invs_noth
[
2
,:,:]
*
np
.
cos
(
rte_invs_noth
[
3
,:,:]
*
np
.
pi
/
180.
)
#blos
rte_data_products
*=
field_stop
[
np
.
newaxis
,
start_row
:
start_row
+
data
.
shape
[
0
],
start_col
:
start_col
+
data
.
shape
[
1
]]
#field stop, set outside to 0
if
out_rte_filename
is
None
:
filename_root
=
str
(
file_path
.
split
(
'.fits'
)[
0
][
-
10
:])
else
:
filename_root
=
out_rte_filename
with
fits
.
open
(
file_path
)
as
hdu_list
:
hdu_list
[
0
].
data
=
rte_data_products
hdu_list
.
writeto
(
out_dir
+
filename_root
+
'_rte_data_products.fits'
,
overwrite
=
True
)
with
fits
.
open
(
file_path
)
as
hdu_list
:
hdu_list
[
0
].
data
=
rte_data_products
[
5
,:,:]
hdu_list
.
writeto
(
out_dir
+
filename_root
+
'_blos_rte.fits'
,
overwrite
=
True
)
with
fits
.
open
(
file_path
)
as
hdu_list
:
hdu_list
[
0
].
data
=
rte_data_products
[
4
,:,:]
hdu_list
.
writeto
(
out_dir
+
filename_root
+
'_vlos_rte.fits'
,
overwrite
=
True
)
with
fits
.
open
(
file_path
)
as
hdu_list
:
hdu_list
[
0
].
data
=
rte_data_products
[
0
,:,:]
hdu_list
.
writeto
(
out_dir
+
filename_root
+
'_Icont_rte.fits'
,
overwrite
=
True
)
printc
(
'--------------------------------------------------------------'
,
bcolors
.
OKGREEN
)
printc
(
f
"------------- CMILOS RTE Run Time:
{
np
.
round
(
time
.
time
()
-
start_time
,
3
)
}
seconds "
,
bcolors
.
OKGREEN
)
printc
(
'--------------------------------------------------------------'
,
bcolors
.
OKGREEN
)
def
pmilos
(
data_f
,
wve_axis_arr
,
data_shape
,
cpos_arr
,
data
,
rte
,
field_stop
,
start_row
,
start_col
,
out_rte_filename
,
out_dir
):
print
(
" "
)
printc
(
'-->>>>>>> RUNNING PMILOS '
,
color
=
bcolors
.
OKGREEN
)
try
:
PMILOS_LOC
=
os
.
path
.
realpath
(
__file__
)
PMILOS_LOC
=
PMILOS_LOC
[:
-
8
]
+
'P-MILOS/'
#11 as hrt_pipe.py is 11 characters -8 if in utils.py
if
os
.
path
.
isfile
(
PMILOS_LOC
+
'pmilos.x'
):
printc
(
"Pmilos executable located at:"
,
PMILOS_LOC
,
color
=
bcolors
.
WARNING
)
else
:
raise
ValueError
(
'Cannot find pmilos:'
,
PMILOS_LOC
)
except
ValueError
as
err
:
printc
(
err
.
args
[
0
],
color
=
bcolors
.
FAIL
)
printc
(
err
.
args
[
1
],
color
=
bcolors
.
FAIL
)
return
wavelength
=
6173.3356
for
scan
in
range
(
int
(
data_shape
[
-
1
])):
start_time
=
time
.
time
()
file_path
=
data_f
[
scan
]
wave_axis
=
wve_axis_arr
[
scan
]
#must invert each scan independently, as cmilos only takes in one dataset at a time
#get wave_axis from the header information of the science scans
if
cpos_arr
[
0
]
==
0
:
shift_w
=
wave_axis
[
3
]
-
wavelength
elif
cpos_arr
[
0
]
==
5
:
shift_w
=
wave_axis
[
2
]
-
wavelength
wave_axis
=
wave_axis
-
shift_w
print
(
'It is assumed the wavelength array is given by the header'
)
#print(wave_axis,color = bcolors.WARNING)
print
(
"Wave axis is: "
,
(
wave_axis
-
wavelength
)
*
1000.
)
print
(
'Saving data into ./P-MILOS/run/data/input_tmp.fits for pmilos RTE input'
)
#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
hdr
=
fits
.
Header
()
primary_hdu
=
fits
.
PrimaryH0DU
(
wave_input
,
header
=
hdr
)
hdul
=
fits
.
HDUList
([
primary_hdu
])
hdul
.
writeto
(
f
'./P-MILOS/run/wavelength_tmp.fits'
,
overwrite
=
True
)
#create input fits file for pmilos
sdata
=
data
[:,:,:,:,
scan
]
hdr
=
fits
.
Header
()
hdr
[
'CTYPE1'
]
=
'HPLT-TAN'
hdr
[
'CTYPE2'
]
=
'HPLN-TAN'
hdr
[
'CTYPE3'
]
=
'STOKES'
#check order of stokes
hdr
[
'CTYPE4'
]
=
'WAVE-GRI'
primary_hdu
=
fits
.
PrimaryHDU
(
sdata
,
header
=
hdr
)
hdul
=
fits
.
HDUList
([
primary_hdu
])
hdul
.
writeto
(
f
'./P-MILOS/run/data/input_tmp.fits'
,
overwrite
=
True
)
#need to change settings for CE or CE+RTE in the pmilos.minit file here
printc
(
f
' ---- >>>>> Inverting data scan number:
{
scan
}
.... '
,
color
=
bcolors
.
OKGREEN
)
cwd
=
os
.
getcwd
()
os
.
chdir
(
"./P-MILOS/run/"
)
cmd
=
"mpiexec -np 16 ../pmilos.x pmilos.minit"
#PMILOS_LOC+"./milos"
if
rte
==
'RTE'
:
cmd
=
"mpiexec -np 16 ../pmilos.x pmilos.minit"
if
rte
==
'CE'
:
cmd
=
"mpiexec -np 16 ../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"
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
:
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
#result has dimensions [rows,cols,13]
"""
PMILOS Output 13 columns
0. eta0 = line-to-continuum absorption coefficient ratio
1. B = magnetic field strength [Gauss]
2. vlos = line-of-sight velocity [km/s]
3. dopp = Doppler width [Angstroms]
4. aa = damping parameter
5. gm = magnetic field inclination [deg]
6. az = magnetic field azimuth [deg]
7. S0 = source function constant
8. S1 = source function gradient
9. mac = macroturbulent velocity [km/s]
10. filling factor of the magnetic component [0-1]
11. Number of iterations performed
12. Chisqr value
"""
rte_data_products
=
np
.
zeros
((
6
,
result
.
shape
[
0
],
result
.
shape
[
1
]))
rte_data_products
[
0
,:,:]
=
result
[:,:,
7
]
+
result
[:,:,
8
]
#continuum
rte_data_products
[
1
,:,:]
=
result
[
1
,:,:]
#b mag strength
rte_data_products
[
2
,:,:]
=
result
[
5
,:,:]
#inclination
rte_data_products
[
3
,:,:]
=
result
[
6
,:,:]
#azimuth
rte_data_products
[
4
,:,:]
=
result
[
2
,:,:]
#vlos
rte_data_products
[
5
,:,:]
=
result
[
1
,:,:]
*
np
.
cos
(
result
[
5
,:,:]
*
np
.
pi
/
180.
)
#blos
rte_data_products
*=
field_stop
[
np
.
newaxis
,
start_row
:
start_row
+
data
.
shape
[
0
],
start_col
:
start_col
+
data
.
shape
[
1
]]
#field stop, set outside to 0
#flipping taken care of for the field stop in the hrt_pipe
if
out_rte_filename
is
None
:
filename_root
=
str
(
file_path
.
split
(
'.fits'
)[
0
][
-
10
:])
else
:
filename_root
=
out_rte_filename
with
fits
.
open
(
file_path
)
as
hdu_list
:
hdu_list
[
0
].
data
=
rte_data_products
hdu_list
.
writeto
(
out_dir
+
filename_root
+
'_rte_data_products.fits'
,
overwrite
=
True
)
with
fits
.
open
(
file_path
)
as
hdu_list
:
hdu_list
[
0
].
data
=
rte_data_products
[
5
,:,:]
hdu_list
.
writeto
(
out_dir
+
filename_root
+
'_blos_rte.fits'
,
overwrite
=
True
)
with
fits
.
open
(
file_path
)
as
hdu_list
:
hdu_list
[
0
].
data
=
rte_data_products
[
4
,:,:]
hdu_list
.
writeto
(
out_dir
+
filename_root
+
'_vlos_rte.fits'
,
overwrite
=
True
)
with
fits
.
open
(
file_path
)
as
hdu_list
:
hdu_list
[
0
].
data
=
rte_data_products
[
0
,:,:]
hdu_list
.
writeto
(
out_dir
+
filename_root
+
'_Icont_rte.fits'
,
overwrite
=
True
)
printc
(
'--------------------------------------------------------------'
,
bcolors
.
OKGREEN
)
printc
(
f
"------------- PMILOS RTE Run Time:
{
np
.
round
(
time
.
time
()
-
start_time
,
3
)
}
seconds "
,
bcolors
.
OKGREEN
)
printc
(
'--------------------------------------------------------------'
,
bcolors
.
OKGREEN
)
\ No newline at end of file
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment