Plot RAD_ALL profiles

# # Reading example of RAD_ALL file.
from astropy.io import fits

gal = 'NGC2916'
t = fits.open(f'{gal}.p_e.RAD_ALL.fits.gz')

### HEADER
# SIMPLE  =                    T / Created with PDL (http://pdl.perl.org)
# BITPIX  =                  -64
# NAXIS   =                    3
# NAXIS1  =                   24
# NAXIS2  =                    6
# NAXIS3  =                   27
# BUNIT   = 'Data Value'
# DESC_23 = 'Lambda  '
# Y_0     = 'R_Re    '
# DESC_2  = 'Sigma_Mass_gas'
# DESC_9  = 'OH_O3N2 '
# DESC_8  = 'OH_t2   '
# DESC_26 = 'NO_N2S2 '
# Y_2     = 'stddev  '
# DESC_12 = 'disp_ssp'
# DESC_0  = 'Sigma_Mass_stars'
# DESC_25 = 'NO_ONS  '
# DESC_5  = 'ZH_MW   '
# DESC_1  = 'Sigma_SFR'
# DESC_13 = 'disp_Ha '
# Y_5     = 'max     '
# Y_3     = 'median  '
# DESC_22 = 'Sigma_Mass_gas_rad'
# DESC_10 = 'Av_gas  '
# DESC_15 = 'vel_ssp '
# DESC_21 = 'Sigma_SFR_ssp'
# DESC_16 = 'SK_ssp  '
# DESC_3  = 'Sigma_Mass_gas_ssp'
# DESC_20 = 'ML      '
# DESC_11 = 'Av_ssp  '
# DESC_19 = 'KIN_Ha  '
# Y_1     = 'mean    '
# DESC_17 = 'SK_Ha   '
# DESC_18 = 'KIN_ssp '
# DESC_7  = 'Age_MW  '
# DESC_4  = 'ZH_LW   '
# DESC_14 = 'vel_Ha  '
# DESC_24 = 'OH_ONS  '
# DESC_6  = 'Age_LW  '
# Y_4     = 'min     '

# t[0].data.shape
# (27, 6, 24)
# 27 parameters, R + 5 statistics (mean, stddev, median, min, max), 24 radial bins
RAD_ALL__psr = t[0].data  # __psr are the dimensions of the vector, i.e., parameters, statistics, radii


# 1st example, the mean (index = 1) Sigma_Mass_stars (index = 0) radial profile is:
Sigma_stars = RAD_ALL__psr[0, 1, :]

# 2nd example, the median (index = 3) stellar metalicity mass weighted (index = 5) radial profile is:
ZH_MW = RAD_ALL__psr[5, 3, :]

# I prefer to work with names instead indices, so you can create you translation dictionary using the header:

h = t[0].header
par_dict = {t[0].header[f'DESC_{i}'].strip(): i for i in range(h['NAXIS1'])}
stats_dict = {t[0].header[f'Y_{i}'].strip(): i for i in range(h['NAXIS2'])}

# par_dict
# {'Sigma_Mass_stars': 0,
#  'Sigma_SFR': 1,
#  'Sigma_Mass_gas': 2,
#  'Sigma_Mass_gas_ssp': 3,
#  'ZH_LW': 4,
#  'ZH_MW': 5,
#  'Age_LW': 6,
#  'Age_MW': 7,
#  'OH_t2': 8,
#  'OH_O3N2': 9,
#  'Av_gas': 10,
#  'Av_ssp': 11,
#  'disp_ssp': 12,
#  'disp_Ha': 13,
#  'vel_Ha': 14,
#  'vel_ssp': 15,
#  'SK_ssp': 16,
#  'SK_Ha': 17,
#  'KIN_ssp': 18,
#  'KIN_Ha': 19,
#  'ML': 20,
#  'Sigma_SFR_ssp': 21,
#  'Sigma_Mass_gas_rad': 22,
#  'Lambda': 23}

# stats_dict
# {'R_Re': 0, 'mean': 1, 'stddev': 2, 'median': 3, 'min': 4, 'max': 5}

# Repeting the above examples using the dictionaries

# 3rd example, the mean Sigma_Mass_stars radial profile:
Sigma_stars = RAD_ALL__psr[:, stats_dict['mean'], par_dict['Sigma_Mass_stars']]

# 4th example, the median ZH_MW radial profile:
ZH_MW = RAD_ALL__psr[:, stats_dict['median'], par_dict['ZH_MW']]

from matplotlib import pyplot as plt

#plotting everything
for _par, _ind in par_dict.items():
    r = RAD_ALL__psr[_ind, stats_dict['R_Re'], :]
    f = plt.figure()
    plt.title(gal)
    plt.xlabel('R [Re]')
    plt.ylabel(_par)
    for _stat in ['mean', 'stddev', 'median', 'min', 'max']:
        plt.plot(r, RAD_ALL__rsp[_ind, stats_dict[_stat], :], label=_stat)
    plt.legend()
    plt.savefig(f'{gal}-{_par}.png')
    plt.close(f)