import math
import numpy as np
from intmodlib import *
import warnings
warnings.filterwarnings("ignore",category=UserWarning)

sigma = 5.670373e-8; #(W m-2 K-4) Stefan-Boltzmann constant
T_sub = 1500.0 #K

########################## Fitter settings #############################

fitter_settings={
'nsteps': 50000,
'nwalkers': 32, #96, #96, #32
'burnin': 20000,
'fit_V_CF': True,
'fit_CP': False, #default: True. If False: closure phases not fitted, if True: closure phases fitted
'fit_DP': False, #default: False
'fit_SED': True, #default: False
#'fit_size': True, #default: False
'deredden_SED': True, #default: False
'A_V': 0.4,
'use_corrflux': False, #default: False. If True: fit correlated fluxes, if False: fit visibilities
'use_chi2': True, #default: True, if False: use likelihood (with parameter log f) for the fitting
'img_size_px': 201, #228
'scale_mas_fixed': 0.3, #0.4, # None, #3.781 #[mas/px] default: None
'zoom_size_mas': None, #15.0, #default: None
'draw_rim': False, #default: True, draw rim for the flared disk model
'pk_version': 2,
'best_fit':'best', #'median' (default), or 'mode', 'best' (best chi2)
'nodata': False, #default: False
'radial_unit': 'au', #'mas' or 'au', default: 'mas'
'color_scheme': 'wl', #color of the symbols on the V vs. baseline plots; None, 'PA', 'wl', default: None 
'img_mode': 'single', #'single' (default), 'avg', or 'pie'
'n_img_avg': 5000, #2000, #number of images to average (or show in a pie image)
'add_img_noise': False, #False (default) or True
'add_radprofile_noise': True, #False (default) or True
'chi2_lim': np.inf, #when making the corner plot, only show the samples with chi2 values better than chi2_lim (default: np.inf)
}

########################### Data settings ##############################

# N band HD 144432
basedir = '/allegro6/matisse/varga/targets3/HD_144432/' #'/Users/jvarga/Data/MATISSE/targets2/HD_144432/' #
#basedir2 = '/allegro6/matisse/varga/test_red/targets_test_spbinN147/HD_144432/'
outputdir = basedir+'/modelling/'

oifits_paths = [ 
#largeAT, good
basedir+'/HD_144432_2019-05-06T07_33_07_Lonly_TARGET_FINALCAL_INT.fits', 
#smallAT, good
basedir+'/HD_144432_2019-09-02T01_19_56_Lonly_TARGET_CHOPPED_FINALCAL_INT.fits', 
#mediumAT, so-so
###basedir+'/HD_144432_2019-09-06T00_31_56_Lonly_TARGET_CHOPPED_FINALCAL_INT.fits', 
#UTs, N band, good
###basedir+'/HD_144432_2020-03-13T08_41_26_N_TARGET_FINALCAL_INT_VISCORRFLUX.fits',
#UTs, N band, so-so
###basedir+'HD_144432_2021-04-03T07_47_34_N_TARGET_FINALCAL_INT_VISCORRFLUX.fits',
#UTs, N band, good
basedir+'HD_144432_2022-03-22T09_11_24_N_TARGET_FINALCAL_INT_VISCORRFLUX.fits',
#UTs, L band, good
basedir+'HD_144432_2022-03-22T09_11_24_Lonly_TARGET_FINALCAL_INT.fits',
#UTs, M band, good
basedir+'HD_144432_2022-03-22T09_11_24_M_TARGET_CHOPPED_FINALCAL_INT.fits',
#medium AT chopped, good
basedir+'HD_144432_2022-03-09T09_10_55_Lonly_TARGET_FINALCAL_INT_0006_AVG.fits',
basedir+'HD_144432_2022-03-09T09_10_55_M_TARGET_FINALCAL_INT_0006_AVG.fits',
#long AT, good
basedir+'HD_144432_2022-04-04T08_23_42_Lonly_TARGET_CHOPPED_FINALCAL_INT.fits',
basedir+'HD_144432_2022-04-04T08_23_42_M_TARGET_CHOPPED_FINALCAL_INT.fits',
#PIONIER, medium AT
basedir+'/2013-06-16_SCI_HD144432_oiDataCalib.fits', 
# large AT
basedir+'/PIONIER.2013-06-07T03 03 55_oidataCalibrated.fits', 
# large AT
#basedir+'/PIONIER.2013-06-10T07 40 34_oidataCalibrated.fits', 
#small AT
basedir+'/PIONIER.2013-07-03T04 09 19_oidataCalibrated.fits', 
basedir+'/GRAVI.2018-03-05T08_29_22.374_singlescivis_singlesciviscalibrated.fits',
#basedir+'/GRAVI.2018-03-05T09_03_37.461_singlescivis_singlesciviscalibrated.fits',
basedir+'/GRAVI.2019-05-24T05_00_42.854_singlescivis_singlesciviscalibrated.fits',
]
'''
oifits_paths = [basedir+'HD_144432_2022-03-22T09_11_24_N_TARGET_FINALCAL_INT_VISCORRFLUX.fits']
'''
oifits_datatypes = [
['int','t3'],
['int','t3'],
['int','t3'],
['int'],#'t3'],
['int'],#'t3'],
['int'],#'t3'],
['int'],#'t3'],
['int'],#'t3'],
['int'],#'t3'],
['int'],#'t3'],
['int'],#'t3'],
['int','t3'],
['int','t3'],
['int','t3','ext1'],
#['int','t3','ext1'],
['int','ext1'],#'t3','ext1'],
] 
'''
oifits_datatypes =[['int','t3']]*10
oifits_datatypes.append(['int','t3','ext1'])
oifits_datatypes.append(['int','t3','ext1'])
oifits_datatypes[2]=['int']
oifits_datatypes[3]=['int']
'''
#oifits_datatypes =[['int']]*1
#oifits_datatypes.append(['int','ext1'])
#oifits_datatypes.append(['int','ext1'])

#print(oifits_datatypes)
#oifits_paths = [basedir+'/HD_144432_2021-04-03T07_47_34_N_TARGET_FINALCAL_INT.fits']
#oifits_datatypes = [['total']]

instrument = 'MATISSE' #MIDI
correct_errors = True
error_terms = {
#'CF_error_percent_lower_limit': 0.08,
#'V_err_lower_limit': 0.03, #not used currently
'V_err_relative_lower_limit': 0.03,
#'V_err_relative_upper_limit': 0.08, #dangerous!!!!!
#'V_err_N_relative_upper_limit': 0.03, #0.003, #give more weigth to the N band data, dangerous!!!!!
'CP_err_lower_limit':  1.0 * math.pi/180.0, #deg
#'CP_err_upper_limit':  5.0 * math.pi/180.0, #deg
#'TF_err_reduce_factor': 0.05,  #!!!!!!!!
#'SED_err_lower_limit_percent': 2.0
}

SEDs = [
    {'path': basedir+'/HD144432_SED_vizier_votable_filtered_cont.vot','format':'votable',
        'includes':['II/297/irc','II/246/out']},#'II/298/fis' #'II/328/allwise', II/246/out 2MASS,'I/305/out','I/350/gaiaedr3'
    {'path': basedir+'/HD144432_MATISSE_LM_SED_with_nans.dat','format':'dat','includes':[]}
] #'/HD144432_MATISSE_LM_SED_with_nans.dat' HD144432_cont_LM_SED_with_nans.dat

stellar_SED_path = basedir+'/stellar_spectrum_PHOENIX_HD_144432.txt' #if you don't want a star,  put ''

misc_spectra = [basedir+'/HD144432_psf.dat'] #,basedir+'/10402662.ss']
sp_labels = ['Spitzer data'] #,'ISOPHOT-S data']

#SED formats: 'votable','dat'
#includes: in case of empty list: [], all data are included

########################## Model settings #############################

#available_models: 
#   'gaussian'
#   'ring_asymm
#   'ring_symm'
#   'temp_grad_asymm' 
#   'temp_grad_symm' 
#   'temp_grad_flared_symm'
#   'TGMdust_symm' 
#   'TGMdust_asymm' 
#   'radimg_symm'
#   'radimg_asymm'
#   'central_1D_symm'
#   'central_1D_asymm'
#   'temp_grad_1D_symm'
#   'temp_grad_1D_asymm'
#   'banana'

model_name = 'TGMdust_symm'
#run_name = 'HKLMN_3zones_symm_withSED_selfconsistenttemp_Mgsilicates_with2Fe_logfract_limited_zone_radii' #_errnoupperlimit' #_except_ri1'
run_name = 'HKLMN_3zones_symm_withSED_tempgrad_Mgsilicates_with2Fe_logfract_logsdens_fix_zone_radii_except_ri1_errnoupperlimit_x3Vchi2weight' #_nocos_i'
#run_name = 'HKLMN_test'
#run_name = 'HKLMN_3zones_symm_withSED_selfconsistenttemp_Mgsilicates_with2Fe_logfract_limited_zone_radii_errnoupperlimit_fitincl'
#run_name = 'asymm_dust_test'
#'HKLMN_gapped_symm_withSED_opt_thin_kappa_after_debug_tempgrad_MgFesilicates_with2Fe' #_4zones_altsil' #_fitall'
#parametrickappa_diffzones'
#run_name = 'HKLN_symm_withSED_opt_thin_kappa_3zones_6species_tg_powerlaw_Nonly' 
#run_name = 'HKLN_gapped_symm_withSED_opt_thin_kappa_crytforstin_pyrmg100out_lt100um_plus_parametrickappa' 
#run_name = 'HKLN_gapped_symm_withSED_opt_thin_kappa_fe_plus_mg100_lt10um'
#run_name = 'HKLN_cont_symm_withSED_opt_thin_kappa_pyr_lt100um_plus_parametrickappa' 
#run_name = 'HKLN_gapped_symm_withSED_opt_thin_kappa_pyr_lt100um_plus_2parametrickappas' 
wl_dic = {
'wavelengths' : [1.7,2.2,3.3,4.75,8.3,11.3,12.3], #[8.3,9.8,11.3,12.5], # [1.7,2.2,3.3,8.3,12.5],#8.3,9.9,12.5], #9.8,11.3 8.5,10.7,12.0], #3.5,
'bandwidths' :  [0.1,0.2,0.2,0.3,0.3, 0.3, 0.3], #[0.3]*4, #[0.1,0.2,0.2,0.4,0.4],#0.2,0.4], #0.2,0.2 ,1.0,1.0,1.0], #(um)
'wl_min' : 1.59, #8.2, #1.59, #1.59, #3.1, #3.1, #8.0, #3.3, #0.2
'wl_max' : 12.4, #12.5, #9.95, #12.5, #12.06, #8.3, #12.0, #12.5, #11.0, #13.0, #10000.0
'wl_mask' : [[1.8,2.1],[2.3,3.2],[3.5,4.7],[4.8,8.27],[8.35,9.79],[9.95,11.26],[11.4,12.2]]
#'wl_max' : 13.0,
#'wl_mask' : [[1.8,2.1],[2.3,3.2],[3.5,4.7],[4.8,8.0]]
#'wl_mask' : [], #[[8.3,8.95],[9.05,9.79],[9.95,10.72],[10.82,11.26],[11.4,12.4]], #[[1.8,2.1],[2.3,3.2],[3.4,8.2],[8.3,9.79],[9.95,12.4]],#[8.3,12.4]], #[3.4,8.2],[8.3,9.79],[9.95,12.4]]
#'wl_mask' : [[1.8,2.1],[2.3,3.2],[3.5,8.2],[8.3,8.95],[9.05,9.79],[9.95,10.72],[10.82,11.26],[11.4,11.99]]
#'wl_mask' : [[1.8,2.1],[2.3,3.2],[3.5,8.44],[8.55,9.97],[10.01,10.95],[11.05,11.99]]
}

# PIONIER data: 1.59-1.79 um

###################### Initial parameter values ########################
#position_angle_init = (76.0)*math.pi/180.0
#axis_ratio_init = 0.87
# stellar parameters HD144432 Vioque+2018
L_star = 9.33 #GAIA DR2 2018: 10.2 #9.33 #L_Sun 
T_star = 7500.0 #7500 #K
R_star = 2.0 #R_Sun #Fairlamb+2015
xc_init = 0.0
yc_init = 0.0
dist_pc = 154.1 #pc #Bailer-Jones GAIA EDR3
position_angle_init = (108.0)*math.pi/180.0 #101.9 #old value: 110.0 #GRAVITY mean: 77.0
axis_ratio_init = 0.79
#position_angle_init = (121.15)*math.pi/180.0 #new MATISSE L-N weighted average of fits
#axis_ratio_init = 0.814 #new!
f_star_init = 0.0 # 0.08
alpha_f_star_init = 1.68
total_flux_star_init = 0.1
#log_f_init = -9.3
total_flux_init = 1.0 #5.5 #Jy, not used in TGMdust model
V_background_init = 0.0
# asymmetry parameters
modulation_amplitude_init = 0.2 #0.51 #0.9
modulation_angle_init = (171.6)*math.pi/180.0 #  109.1
r_mod_max_init = 10.0 #0.6/dist_pc*1000.0 #[mas]
sigma_r_mod_init = 14.0 #0.2/dist_pc*1000.0 #[mas]

#PIONIER hlr = 1.5 mas = 0.23 au

# disk parameters
r_taper_au_init = 30.0 #30.0 #30.0 #should be 40.0 #au #MATISSE pinhole diameter, 2 lambda/D, on UTs at 10 um it is 0.5'' = 78 au
Ti1_init = 1729.2  #1578.1
#withC init
ri1_init = 0.17 #0.16 #0.23 #0.17 #0.24 #0.17 #0.21 #0.17 #au
ri2_init = 1.26 #1.21 #1.29 #0.87 #0.49 #1.33 #0.4 #1.19 #0.55 #au
ri3_init = 4.13 ######4.49 #3.95 #4.11 #4.27 #3.99
p1_init = -0.788 #-1.073 #-8.1 #-13.33 #8.486#-0.755 #0.152 #-0.44 #-8.12 #-0.44 #-2.325 #-3.413 #
p2_init = -1.166 #########0.488 # 0.717 ##0.681 #0.134 #0.074 #0.676 #0.074 #0.0 #2.5 #-2.81 #-0.372 #-1.212 #
p3_init = -62.256 #-13.14 #-4.8 #-1.83 # #-8.921 #-8.146 #-8.441 #1.0
#dr_smooth_init = 0.07 
#C_i_init = 0.0 #0.0008 #0.00008 #0.0005
#C_o_init = 0.0 #0.00008 #0.08 #0.00002
tg_sur_init = -0.529 #-0.748 #-0.614 #-0.575 #-0.573 #-0.596 #-0.573

inner_rim_boost_factor_init = 1.0
 
kappa0_init = 1100.0
lambda_break_init = 4.18 #2.13 #0.9 #2.13 #um
p0_kappa_init = 1.991#0.2 #-0.04
p1_kappa_init = 2.394 #0.075 #2.1

N_zones = 3
norm1o = -3.925 #np.log10(0.000134) #0.000168 #0.00048 #0.0015998  #0.000539
norm2o = -3.469 #np.log10(0.000277) #0.000246#3.1849267e-05 #0.000123 #0.000043
norm3o = -0.032 #np.log10(1.904) ##########0.0634 # 0.001754 #0.001496 #0.001249

'''
#withFe init
ri1_init = 0.21 #0.24 #0.17 #0.21 #0.17 #au
ri2_init = 1.0 #0.49 #1.33 #0.4 #1.19 #0.55 #au
ri3_init = 4.3
#ri4_init = 4.11 #4.27 #3.99
p1_init = -2.11 #-0.755 #0.152 #-0.44 #-8.12 #-0.44 #-2.325 #-3.413 #
p2_init = 0.218
p3_init = -13.36 #0.681 #0.134 #0.074 #0.676 #0.074 #0.0 #2.5 #-2.81 #-0.372 #-1.212 #
#p4_init = -8.921 #-8.146 #-8.441 #1.0
#dr_smooth_init = 0.07 
#C_i_init = 0.0 #0.0008 #0.00008 #0.0005
#C_o_init = 0.0 #0.00008 #0.08 #0.00002
tg_sur_init = -0.538 #-0.575 #-0.573 #-0.596 #-0.573

kappa0_init = 1100.0
lambda_break_init = 4.18 #2.13 #0.9 #2.13 #um
p0_kappa_init = 1.991#0.2 #-0.04
p1_kappa_init = 2.394 #0.075 #2.1

N_zones = 3
norm1o = 0.000763 #0.00021  #0.000531 #0.000194  #0.000434495 #0.0084088 #0.00038523 #0.003326 #0.000424
norm2o = 0.000278 #0.0001
norm3o = 0.01639 #0.000159 #0.000044 #0.000343 #0.001065 #0.00029114 #0.000055 #0.000263
#norm4o = 0.007233 #0.000838 #0.008546 #0.006859 # 0.00014667 #0.00068 #0.00013
'''
norm1 = 1.0
norm2 = 1.0
norm3 = 1.0
# mineralogy parameters

#with C
'''
sdens_dic={
'sdens_z1_1_init' : norm1*0.18, # [g cm-2] 
'sdens_z1_2_init' : norm1*0.05, 
'sdens_z1_3_init' : norm1*0.10, 
'sdens_z1_4_init' : norm1*0.02, 
'sdens_z1_5_init' : norm1*0.51, 
'sdens_z1_6_init' : norm1*0.14, 
'sdens_z2_1_init' : norm2*0.49, 
'sdens_z2_2_init' : norm2*0.22, 
'sdens_z2_3_init' : norm2*0.00, 
'sdens_z2_4_init' : norm2*0.01, 
'sdens_z2_5_init' : norm2*0.13, 
'sdens_z2_6_init' : norm2*0.14, 
'sdens_z3_1_init' : norm3*0.16,
'sdens_z3_2_init' : norm3*0.04, 
'sdens_z3_3_init' : norm3*0.00,
'sdens_z3_4_init' : norm3*0.01, 
'sdens_z3_5_init' : norm3*0.05, 
'sdens_z3_6_init' : norm3*0.74, 
}
'''
'''
#with Fe new fit
sdens_dic={
'sdens_z1_1_init' : norm1*0.01, # [g cm-2] 
'sdens_z1_2_init' : norm1*0.00, 
'sdens_z1_3_init' : norm1*0.02, 
'sdens_z1_4_init' : norm1*0.04, 
'sdens_z1_5_init' : norm1*0.01, 
'sdens_z1_6_init' : norm1*0.00,#0.416,#0.43,
'sdens_z1_7_init' : norm1*0.01,
'sdens_z1_8_init' : norm1*0.00,
'sdens_z1_9_init' : norm1*0.14,
'sdens_z1_10_init' : norm1*0.32,
'sdens_z1_11_init' : norm1*0.47,
'sdens_z2_1_init' : norm2*0.00, # [g cm-2] 
'sdens_z2_2_init' : norm2*0.14, 
'sdens_z2_3_init' : norm2*0.07, 
'sdens_z2_4_init' : norm2*0.00, 
'sdens_z2_5_init' : norm2*0.0, 
'sdens_z2_6_init' : norm2*0.0,#0.416,#0.43,
'sdens_z2_7_init' : norm2*0.0,
'sdens_z2_8_init' : norm2*0.0,
'sdens_z2_9_init' : norm2*0.04,
'sdens_z2_10_init' : norm2*0.04,
'sdens_z2_11_init' : norm2*0.71,
'sdens_z3_1_init' : norm3*0.0, 
'sdens_z3_2_init' : norm3*0.0, 
'sdens_z3_3_init' : norm3*0.0, 
'sdens_z3_4_init' : norm3*0.0, 
'sdens_z3_5_init' : norm3*0.0, 
'sdens_z3_6_init' : norm3*0.0,#0.66,# 0.55, 
'sdens_z3_7_init' : norm3*0.0,
'sdens_z3_8_init' : norm3*0.0,
'sdens_z3_9_init' : norm3*0.0,
'sdens_z3_10_init' : norm3*0.01,
'sdens_z3_11_init' : norm3*0.99,
}
'''
'''
#with C new fit
sdens_dic={
'sdens_z1_1_init' : norm1*0.00, # [g cm-2] 
'sdens_z1_2_init' : norm1*0.28, 
'sdens_z1_3_init' : norm1*0.03, 
'sdens_z1_4_init' : norm1*0.09, 
'sdens_z1_5_init' : norm1*0.03, 
'sdens_z1_6_init' : norm1*0.00,#0.416,#0.43,
'sdens_z1_7_init' : norm1*0.01,
'sdens_z1_8_init' : norm1*0.00,
'sdens_z1_9_init' : norm1*0.49,
'sdens_z1_10_init' : norm1*0.06,
'sdens_z1_11_init' : norm1*0.01,
'sdens_z2_1_init' : norm2*0.00, # [g cm-2] 
'sdens_z2_2_init' : norm2*0.18, 
'sdens_z2_3_init' : norm2*0.44, 
'sdens_z2_4_init' : norm2*0.07, 
'sdens_z2_5_init' : norm2*0.02, 
'sdens_z2_6_init' : norm2*0.07,#0.416,#0.43,
'sdens_z2_7_init' : norm2*0.00,
'sdens_z2_8_init' : norm2*0.07,
'sdens_z2_9_init' : norm2*0.00,
'sdens_z2_10_init' : norm2*0.15,
'sdens_z2_11_init' : norm2*0.01,
'sdens_z3_1_init' : norm3*0.0, 
'sdens_z3_2_init' : norm3*0.0, 
'sdens_z3_3_init' : norm3*0.46, 
'sdens_z3_4_init' : norm3*0.0, 
'sdens_z3_5_init' : norm3*0.01, 
'sdens_z3_6_init' : norm3*0.08,#0.66,# 0.55, 
'sdens_z3_7_init' : norm3*0.0,
'sdens_z3_8_init' : norm3*0.06,
'sdens_z3_9_init' : norm3*0.0,
'sdens_z3_10_init' : norm3*0.38,
'sdens_z3_11_init' : norm3*0.01,
}
'''
'''
#with Mg silicates, with Fe new fit 2023 May
sdens_dic={
'sdens_z1_1_init' : norm1*np.log10(0.00000000001), # [g cm-2] 0.18
'sdens_z1_2_init' : norm1*np.log10(0.0599), 
'sdens_z1_3_init' : norm1*np.log10(0.00000000001), 
'sdens_z1_4_init' : norm1*np.log10(0.0305), 
'sdens_z1_5_init' : norm1*np.log10(0.5109),
'sdens_z1_6_init' : norm1*np.log10(0.39/2.0), #(-1.337),
'sdens_z1_7_init' : norm1*np.log10(0.39/2.0), #0.542,
'sdens_z2_1_init' : norm2*np.log10(0.261),
'sdens_z2_2_init' : norm2*np.log10(0.198), 
'sdens_z2_3_init' : norm2*np.log10(0.00000000001), 
'sdens_z2_4_init' : norm2*np.log10(0.0129), 
'sdens_z2_5_init' : norm2*np.log10(0.168), 
'sdens_z2_6_init' : norm2*np.log10(0.36/2.0), #(-2.385),
'sdens_z2_7_init' : norm2*np.log10(0.36/2.0), #(0.373),
'sdens_z3_1_init' : norm3*np.log10(0.0449), 
'sdens_z3_2_init' : norm3*np.log10(0.0235), 
'sdens_z3_3_init' : norm3*np.log10(0.00000000001), 
'sdens_z3_4_init' : norm3*np.log10(0.00492), 
'sdens_z3_5_init' : norm3*np.log10(0.00644), 
'sdens_z3_6_init' : norm3*np.log10(0.92/2.0), #(-3.29),
'sdens_z3_7_init' : norm3*np.log10(0.92/2.0), #(-1.39),
}
'''

#with Mg silicates, with Fe new fit 2023 before May, I use this in the paper
sdens_dic={
'sdens_z1_1_init' : norm1*np.log10(0.18), # [g cm-2] 0.18
'sdens_z1_2_init' : norm1*np.log10(0.0363), 
'sdens_z1_3_init' : norm1*np.log10(0.0596), 
'sdens_z1_4_init' : norm1*np.log10(0.0146), 
'sdens_z1_5_init' : norm1*np.log10(0.425),
'sdens_z1_6_init' : norm1*(-2.496),# (-0.681),
'sdens_z1_7_init' : norm1*( 0.916),#(-1.481),
'sdens_z2_1_init' : norm2*np.log10(0.211),
'sdens_z2_2_init' : norm2*np.log10(0.168), 
'sdens_z2_3_init' : norm2*np.log10(0.0000001), 
'sdens_z2_4_init' : norm2*np.log10(0.0036), 
'sdens_z2_5_init' : norm2*np.log10(0.0956), 
'sdens_z2_6_init' : norm2*(-0.306),#(-1.280),
'sdens_z2_7_init' : norm2*(-0.940),#(-2.664),
'sdens_z3_1_init' : norm3*np.log10(0.0433), 
'sdens_z3_2_init' : norm3*np.log10(0.0000001), 
'sdens_z3_3_init' : norm3*np.log10(0.0000001), 
'sdens_z3_4_init' : norm3*np.log10(0.0000001), 
'sdens_z3_5_init' : norm3*np.log10(0.010), 
'sdens_z3_6_init' : norm3*(-2.796),#(-3.904),
'sdens_z3_7_init' : norm3*(-5.268),#(-4.985),
}

'''
#with Mg silicates, with C new fit 2023 before May, I use this in the paper
sdens_dic={
'sdens_z1_1_init' : norm1*np.log10(0.18), # [g cm-2] 0.18
'sdens_z1_2_init' : norm1*np.log10(0.0363), 
'sdens_z1_3_init' : norm1*np.log10(0.0596), 
'sdens_z1_4_init' : norm1*np.log10(0.0146), 
'sdens_z1_5_init' : norm1*np.log10(0.425),
'sdens_z1_6_init' : norm1*(-0.681),
'sdens_z1_7_init' : norm1*(-1.481),
'sdens_z2_1_init' : norm2*np.log10(0.211),
'sdens_z2_2_init' : norm2*np.log10(0.168), 
'sdens_z2_3_init' : norm2*np.log10(0.0000001), 
'sdens_z2_4_init' : norm2*np.log10(0.0036), 
'sdens_z2_5_init' : norm2*np.log10(0.0956), 
'sdens_z2_6_init' : norm2*(-1.280),
'sdens_z2_7_init' : norm2*(-2.664),
'sdens_z3_1_init' : norm3*np.log10(0.0433), 
'sdens_z3_2_init' : norm3*np.log10(0.0000001), 
'sdens_z3_3_init' : norm3*np.log10(0.0000001), 
'sdens_z3_4_init' : norm3*np.log10(0.0000001), 
'sdens_z3_5_init' : norm3*np.log10(0.010), 
'sdens_z3_6_init' : norm3*(-3.904),
'sdens_z3_7_init' : norm3*(-4.985),
}
'''
'''
#with MgFe silicates, with Fe new fit
sdens_dic={
'sdens_z1_1_init' : norm1*0.18, # [g cm-2] 
'sdens_z1_2_init' : norm1*0.0, 
'sdens_z1_3_init' : norm1*0.0, 
'sdens_z1_4_init' : norm1*0.0, 
'sdens_z1_5_init' : norm1*0.01, 
'sdens_z1_6_init' : norm1*0.0,#0.416,#0.43,
'sdens_z1_7_init' : norm1*0.0,
'sdens_z1_8_init' : norm1*0.02,
'sdens_z1_9_init' : norm1*0.78,
'sdens_z1_10_init' : norm1*0.01,
#'sdens_z1_11_init' : norm1*0.01,
'sdens_z2_1_init' : norm2*0.14, # [g cm-2] 
'sdens_z2_2_init' : norm2*0.0, 
'sdens_z2_3_init' : norm2*0.29, 
'sdens_z2_4_init' : norm2*0.0, 
'sdens_z2_5_init' : norm2*0.0, 
'sdens_z2_6_init' : norm2*0.0,#0.416,#0.43,
'sdens_z2_7_init' : norm2*0.0,
'sdens_z2_8_init' : norm2*0.0,
'sdens_z2_9_init' : norm2*0.56,
'sdens_z2_10_init' : norm2*0.01,
#'sdens_z2_11_init' : norm2*0.01,
'sdens_z3_1_init' : norm3*0.0, 
'sdens_z3_2_init' : norm3*0.0, 
'sdens_z3_3_init' : norm3*0.07, 
'sdens_z3_4_init' : norm3*0.0, 
'sdens_z3_5_init' : norm3*0.0, 
'sdens_z3_6_init' : norm3*0.0,#0.66,# 0.55, 
'sdens_z3_7_init' : norm3*0.0,
'sdens_z3_8_init' : norm3*0.0,
'sdens_z3_9_init' : norm3*0.92,
'sdens_z3_10_init' : norm3*0.01,
#'sdens_z3_11_init' : norm3*0.01,
}
'''
'''
#with MgFe silicates, with C new fit
sdens_dic={
'sdens_z1_1_init' : norm1*0.66, # [g cm-2] 
'sdens_z1_2_init' : norm1*0.0, 
'sdens_z1_3_init' : norm1*0.0, 
'sdens_z1_4_init' : norm1*0.0, 
'sdens_z1_5_init' : norm1*0.07, 
'sdens_z1_6_init' : norm1*0.0,#0.416,#0.43,
'sdens_z1_7_init' : norm1*0.0,
'sdens_z1_8_init' : norm1*0.07,
'sdens_z1_9_init' : norm1*0.18,
'sdens_z1_10_init' : norm1*0.01,
#'sdens_z1_11_init' : norm1*0.01,
'sdens_z2_1_init' : norm2*0.37, # [g cm-2] 
'sdens_z2_2_init' : norm2*0.0, 
'sdens_z2_3_init' : norm2*0.56, 
'sdens_z2_4_init' : norm2*0.0, 
'sdens_z2_5_init' : norm2*0.0, 
'sdens_z2_6_init' : norm2*0.0,#0.416,#0.43,
'sdens_z2_7_init' : norm2*0.0,
'sdens_z2_8_init' : norm2*0.0,
'sdens_z2_9_init' : norm2*0.05,
'sdens_z2_10_init' : norm2*0.01,
#'sdens_z2_11_init' : norm2*0.01,
'sdens_z3_1_init' : norm3*0.0, 
'sdens_z3_2_init' : norm3*0.0, 
'sdens_z3_3_init' : norm3*0.60, 
'sdens_z3_4_init' : norm3*0.0, 
'sdens_z3_5_init' : norm3*0.0, 
'sdens_z3_6_init' : norm3*0.0,#0.66,# 0.55, 
'sdens_z3_7_init' : norm3*0.0,
'sdens_z3_8_init' : norm3*0.0,
'sdens_z3_9_init' : norm3*0.39,
'sdens_z3_10_init' : norm3*0.01,
#'sdens_z3_11_init' : norm3*0.01,
}
'''

############################# Model setup ##############################

# - define models with parameters
if model_name.startswith('temp_grad_flared'):
    pa_range = [0.0,2.0*np.pi]
else:
    pa_range = [0.0,1.0*np.pi]

F = False
T = True
FI = False
FO = False
TT = False
TU = True
TR = False
# arguments of modelparameter init function: param_key,fitted,priors,init_value,param_range,label,unit,print_format,print_factor
# geometry and flux parameters
parameter_dic = {
'common': [
modelparameter('position_angle',      F ,[0.0*np.pi/180.0,180.0*np.pi/180.0],position_angle_init,None,r'\theta',r'\degree',r'\mathrm{rad}','%.1f',180.0/np.pi),
modelparameter('axis_ratio',          F ,[0.3,1.0],axis_ratio_init,None,r'\cos\,i','',None,'%.2f',1.0),
modelparameter('f_star',              F ,[0.0,1.0],f_star_init,[f_star_init*0.95,f_star_init*1.05],'f_\mathrm{star}','',None,'%.2f',1.0),
modelparameter('alpha_f_star',        F ,[0.0,3.0],alpha_f_star_init,None,r'\alpha_{f\mathrm{star}}','',None,'%.2f',1.0),
modelparameter('total_flux_star',     F ,[0.01,100.0],total_flux_star_init,None,r'F_\mathrm{tot,star}','\mathrm{Jy}',None,'%.2f',1.0),
#modelparameter('log_f',              F ,[-10.0,2.0],log_f_init,[log_f_init-1.0,log_f_init+1.0],r'\log\,f','',None,'%.2f',1.0),
modelparameter('total_flux',          F ,[0.01,150.0],total_flux_init,None,r'F_\mathrm{tot}','\mathrm{Jy}',None,'%.1f',1.0),
modelparameter('V_background',        F ,[0.0,1.0],V_background_init,None,r'V_\mathrm{bg}','',None,'%.2f',1.0),
modelparameter('xc',                  F ,[-7.0,7.0],xc_init,[xc_init-1.0,xc_init+1.0],r'x_c','\mathrm{mas}',None,'%.2f',1.0),
modelparameter('yc',                  F ,[-7.0,7.0],yc_init,[yc_init-1.0,yc_init+1.0],r'y_c','\mathrm{mas}',None,'%.2f',1.0),
],
# parameters of the asymmetry
'common_asymm': [
modelparameter('modulation_amplitude',F,[0.0,1.0],modulation_amplitude_init,None,r'A_\mathrm{mod}','',None,'%.2f',1.0),
modelparameter('modulation_angle',    F,[0.0,2.0*np.pi],modulation_angle_init,None,r'\phi_\mathrm{mod}',r'\degree',r'\mathrm{rad}','%.1f',180.0/np.pi),
modelparameter('modulation_radius',   F,[0.01,100.0],r_mod_max_init,None,r'r_\mathrm{mod,max}','\mathrm{mas}',None,'%.2f',1.0),
modelparameter('modulation_radial_width',F ,[0.01,60.0],sigma_r_mod_init,None,r'\sigma_\mathrm{r,mod}','\mathrm{mas}',None,'%.2f',1.0)
],
#'radimg': [
'TGMdust': [
# stellar parameters
modelparameter('L_star',              F ,[0.1,1000.0],L_star,None,r'L_*',r'\mathrm{L}_\odot',None,'%.1f',1.0),
modelparameter('T_star',              F ,[1000.0,25000.0],T_star,None,r'T_*',r'\mathrm{K}',None,'%.1f',1.0),
modelparameter('R_star',              F ,[0.1,1000.0],R_star,None,r'R_*',r'\mathrm{R}_\odot',None,'%.1f',1.0),
modelparameter('dist_pc',             F ,[1.0,2000.0],dist_pc,None,r'd',r'\mathrm{pc}',None,'%.1f',1.0),
# disk parameters
modelparameter('r_taper_au',          F ,[1.0,1000.0],r_taper_au_init,None,r'R_\mathrm{taper}','\mathrm{au}',None,'%.2f',1.0),
modelparameter('T_i1',                T ,[1200.0,2000.0],Ti1_init,None,r'T_\mathrm{in}','\mathrm{K}',None,'%.1f',1.0),
modelparameter('r_i1',                T ,[0.15,0.25],ri1_init,None,r'R_1','\mathrm{au}',None,'%.2f',1.0),
modelparameter('r_i2',                F ,[1.04-0.023,1.04+0.058],ri2_init,None,r'R_2','\mathrm{au}',None,'%.2f',1.0),
modelparameter('r_i3',                F ,[4.172-0.004,4.172+0.843],ri3_init,None,r'R_3','\mathrm{au}',None,'%.2f',1.0),
#modelparameter('r_i4',                T ,[2.1,8.0],ri4_init,None,r'R_\mathrm{i4}','\mathrm{au}',None,'%.2f',1.0),
#modelparameter('r_gap_in',            F ,[0.01,10.0],r_gap_in_init,None,r'R_\mathrm{gap,in}','\mathrm{au}',None,'%.2f',1.0),
#modelparameter('dr_smooth',           T ,[0.0,1.0],dr_smooth_init,None,r'dr_\mathrm{smooth}','\mathrm{au}',None,'%.2f',1.0),
modelparameter('tg_sur',              T ,[-1.5,-0.3],tg_sur_init,None,r'q','',None,'%.2f',1.0), #temperature gradient power law
modelparameter('p1',                  T ,[-10.0,3.0],p1_init,None,r'p_1','',None,'%.2f',1.0), 
modelparameter('p2',                  T ,[-5.0,5.0],p2_init,None,r'p_2','',None,'%.2f',1.0), 
modelparameter('p3',                  T ,[-100.0,4.0],p3_init,None,r'p_3','',None,'%.2f',1.0), 
#modelparameter('p4',                  T ,[-16.0,4.0],p4_init,None,r'p_4','',None,'%.3f',1.0), 
modelparameter('n1',                  T ,[-9.0,1.0],norm1o,None,r'\log_{10} \Sigma_{0,1}',r'\mathrm{g\ cm}^{-2}',None,'%.6f',1.0), 
modelparameter('n2',                  T ,[-9.0,1.0],norm2o,None,r'\log_{10} \Sigma_{0,2}',r'\mathrm{g\ cm}^{-2}',None,'%.6f',1.0), 
modelparameter('n3',                  T ,[-9.0,2.0],norm3o,None,r'\log_{10} \Sigma_{0,3}',r'\mathrm{g\ cm}^{-2}',None,'%.6f',1.0), 
#modelparameter('inner_rim_boost_factor',T ,[0.1,20.0],inner_rim_boost_factor_init,None,r'f_\mathrm{boost,rim}','',None,'%.2f',1.0), 
#modelparameter('n4',                  T ,[0.0,0.1],norm4o,None,r'\Sigma_4',r'\mathrm{g\ cm}^{-2}',None,'%.6f',1.0), 
#modelparameter('r_break',             T ,[0.01,100.0],r_break_init,None,r'r_\mathrm{break}',r'\mathrm{au}',None,'%.2f',1.0),
#modelparameter('r_qjump',             F ,[0.01,100.0],r_qjump_init,None,r'r_\mathrm{qjump}',r'\mathrm{au}',None,'%.2f',1.0),
# modelparameter('kappa0',              F, [0.0,10000.0],kappa0_init,None,r'\kappa_0',r'\mathrm{cm}^{2}\ \mathrm{g}^{-1}',None,'%.3f',1.0),
# modelparameter('lambda_break_z1',        T, [0.01,100.0],lambda_break_init,None,r'\lambda_{\kappa_\mathrm{break,z1}}',r'\mu\mathrm{m}',None,'%.2f',1.0),
# modelparameter('p0_kappa_z1',            T, [-3.0,4.0],p0_kappa_init,None,r'p_{0,\kappa,z1}','',None,'%.3f',1.0),
# modelparameter('p1_kappa_z1',            T, [-3.0,4.0],p1_kappa_init,None,r'p_{1,\kappa,z1}','',None,'%.3f',1.0),
# modelparameter('lambda_break_z2',        T, [0.01,100.0],lambda_break_init,None,r'\lambda_{\kappa_\mathrm{break,z2}}',r'\mu\mathrm{m}',None,'%.2f',1.0),
# modelparameter('p0_kappa_z2',            T, [-3.0,4.0],p0_kappa_init,None,r'p_{0,\kappa,z2}','',None,'%.3f',1.0),
# modelparameter('p1_kappa_z2',            T, [-3.0,4.0],p1_kappa_init,None,r'p_{1,\kappa,z2}','',None,'%.3f',1.0),
# modelparameter('lambda_break_z3',        T, [0.01,100.0],lambda_break_init,None,r'\lambda_{\kappa_\mathrm{break,z3}}',r'\mu\mathrm{m}',None,'%.2f',1.0),
# modelparameter('p0_kappa_z3',            T, [-3.0,4.0],p0_kappa_init,None,r'p_{0,\kappa,z3}','',None,'%.3f',1.0),
# modelparameter('p1_kappa_z3',            T, [-3.0,4.0],p1_kappa_init,None,r'p_{1,\kappa,z3}','',None,'%.3f',1.0),
#modelparameter('okappa0',              F, [0.0,10000.0],kappa0_init,None,r'\kappa_{o0}',r'\mathrm{cm}^{2}\ \mathrm{g}^{-1}',None,'%.3f',1.0),
#modelparameter('olambda_break',        T, [0.01,10000.0],lambda_break_init,None,r'\lambda_{\kappa_\mathrm{o,break}}',r'\mu\mathrm{m}',None,'%.2f',1.0),
#modelparameter('op0_kappa',            T, [-3.0,3.0],p0_kappa_init,None,r'p_{o0,\kappa}','',None,'%.3f',1.0),
#modelparameter('op1_kappa',            T, [-3.0,3.0],p1_kappa_init,None,r'p_{o1,\kappa}','',None,'%.3f',1.0),
]
}

ks = []
'''
ks.append('dustkappa_pure_Fe.dat')
ks.append('dustkappa_default.dat' )
'''

#ks.append('Q_Am_Mgolivine_Jae_DHS_f0.7_rv0.1.dat')
ks.append('Q_Am_Mgolivine_Jae_DHS_f0.7_rv2.0.dat')
ks.append('Q_Am_Mgpyroxene_Dor_DHS_f0.7_rv0.1.dat')
ks.append('Q_Am_Mgpyroxene_Dor_DHS_f0.7_rv2.0.dat')
ks.append('Q_Fo_Sogawa_DHS_f1.0_rv0.1.dat')
#ks.append('Q_Fo_Sogawa_DHS_f1.0_rv2.0.dat')
#ks.append('Q_En_Jaeger_DHS_f1.0_rv0.1.dat')
#ks.append('Q_En_Jaeger_DHS_f1.0_rv2.0.dat')
ks.append('Q_En_Jaeger_DHS_f1.0_rv5.0.dat')
ks.append('Q_iron_0.10um_dhs_0.70.dat')
ks.append('Q_iron_2.00um_dhs_0.70.dat')
#ks.append('Q_c-z_0.10um_dhs_0.70.dat')
#ks.append('Q_c-z_2.00um_dhs_0.70.dat')
'''
ks.append('Q_ol-mg40_0.10um_dhs_0.70.dat')
ks.append('Q_ol-mg40_2.00um_dhs_0.70.dat')
ks.append('Q_pyr-mg40_0.10um_dhs_0.70.dat')
ks.append('Q_pyr-mg40_2.00um_dhs_0.70.dat')
ks.append('Q_fay_0.10um_dhs_0.99.dat')
ks.append('Q_fay_2.00um_dhs_0.99.dat')
ks.append('Q_En_Jaeger_DHS_f1.0_rv0.1.dat')
ks.append('Q_En_Jaeger_DHS_f1.0_rv2.0.dat')
ks.append('Q_iron_0.10um_dhs_0.70.dat')
ks.append('Q_iron_2.00um_dhs_0.70.dat')
#ks.append('Q_c-z_0.10um_dhs_0.70.dat')
#ks.append('Q_c-z_2.00um_dhs_0.70.dat')
'''
'''
ks.append('Q_Am_Mgolivine_Jae_DHS_f0.7_rv2.0.dat')
ks.append('Q_Am_Mgpyroxene_Dor_DHS_f0.7_rv0.1.dat')
ks.append('Q_Am_Mgpyroxene_Dor_DHS_f0.7_rv2.0.dat')
ks.append('Q_Fo_Sogawa_DHS_f1.0_rv0.1.dat')
ks.append('Q_En_Jaeger_DHS_f1.0_rv5.0.dat')
ks.append('Q_fe-c_0.10-0.30um.dat')
#ks.append('Q_amorph_c_rv0.1.dat')
#ks.append('parametric')
'''
N_dust = len(ks)

parkappa_idx = None
for i in range(N_dust):
    if ks[i] == 'parametric':
        parkappa_idx = i
#print(parkappa_idx)

for i in range(N_zones):
    for j in range(N_dust):
        #parameter_dic['TGMdust'].append(modelparameter('sdens_z%d_%d'%(i+1,j+1),    \         
        #    F, [0.0,10.0],sdens_dic['sdens_z%d_%d_init'%(i+1,j+1)],None,r'\Sigma_\mathrm{z%d\ %d}'%(i+1,j+1),\
        #    r'\mathrm{g\ cm}^{-2}',None,'%.3f',1.0))
        if j == parkappa_idx:
            fit = T
        else:
            fit = F
        # parameter_dic['TGMdust'].append(modelparameter('sdens_z%d_%d'%(i+1,j+1),             
            # fit, [0.0,1.0],sdens_dic['sdens_z%d_%d_init'%(i+1,j+1)],None,r'f_\mathrm{z%d\ %d}'%(i+1,j+1),\
            # '',None,'%.3f',1.0))
        init_val = sdens_dic['sdens_z%d_%d_init'%(i+1,j+1)]
        if 'iron' in ks[j] or 'c-z' in ks[j]:
            range_min = init_val-3.0
            range_max = init_val+2.0
            if range_max > 1.5: #15.0:
                range_max = 1.5 #15.0
            text=''
            if j+1 == 6:
                text = 'small'
            if j+1 == 7:
                text = 'large'
            parameter_dic['TGMdust'].append(modelparameter('sdens_z%d_%d'%(i+1,j+1),             
                T, [range_min,range_max],init_val,None,r'c_\mathrm{%d\ %s}'%(i+1,text),\
                '',None,'%.3f',1.0))
        else:
            range_min = 0.75*init_val
            range_max = 1.25*init_val
            parameter_dic['TGMdust'].append(modelparameter('sdens_z%d_%d'%(i+1,j+1),             
                F, [range_min,range_max],init_val,None,r'c_\mathrm{z%d\ %d}'%(i+1,j+1),\
                '',None,'%.3f',1.0))
        print('sdens_z%d_%d'%(i+1,j+1),init_val,range_min,range_max)

#modelparameter('sdens_z1_2',             T, [0.0,100.0],sdens_z1_2_init,None,r'\Sigma_\mathrm{in2}',r'\mathrm{g\ cm}^{-2}',None,'%.6f',1.0),
# modelparameter('sdens_z1_3',             T, [0.0,100.0],sdens_z1_3_init,None,r'\Sigma_\mathrm{in3}',r'\mathrm{g\ cm}^{-2}',None,'%.6f',1.0),
#modelparameter('sdens_z2_1',             T, [0.0,10.0],sdens_z2_1_init,None,r'\Sigma_\mathrm{m1}',r'\mathrm{g\ cm}^{-2}',None,'%.6f',1.0),
# modelparameter('sdens_z2_2',             T, [0.0,100.0],sdens_z2_2_init,None,r'\Sigma_\mathrm{m2}',r'\mathrm{g\ cm}^{-2}',None,'%.6f',1.0),
# modelparameter('sdens_z2_3',             T, [0.0,100.0],sdens_z2_3_init,None,r'\Sigma_\mathrm{m3}',r'\mathrm{g\ cm}^{-2}',None,'%.6f',1.0),
#modelparameter('sdens_z3_1',             T, [0.0,10.0],sdens_z3_1_init,None,r'\Sigma_\mathrm{out1}',r'\mathrm{g\ cm}^{-2}',None,'%.6f',1.0),
# modelparameter('sdens_z3_2',             T, [0.0,100.0],sdens_z3_2_init,None,r'\Sigma_\mathrm{out2}',r'\mathrm{g\ cm}^{-2}',None,'%.6f',1.0),
# modelparameter('sdens_z3_3',             T, [0.0,100.0],sdens_z3_3_init,None,r'\Sigma_\mathrm{out3}',r'\mathrm{g\ cm}^{-2}',None,'%.6f',1.0),


T=True
# q1= 'Q_Am_Mgolivine_Jae_DHS_f0.7_rv0.1.dat'
# q2 = 'Q_Am_Mgolivine_Jae_DHS_f0.7_rv2.0.dat'
# q3 = 'Q_Am_Mgpyroxene_Dor_DHS_f0.7_rv0.1.dat'
# q4 = 'Q_Am_Mgpyroxene_Dor_DHS_f0.7_rv2.0.dat'
# q5 = 'Q_Fo_Sogawa_DHS_f1.0_rv0.1.dat'
# q6 = 'Q_Fo_Sogawa_DHS_f1.0_rv2.0.dat'
# q7 = 'Q_En_Jaeger_DHS_f1.0_rv0.1.dat'
# q8 = 'Q_En_Jaeger_DHS_f1.0_rv2.0.dat'
# q9 = 'Q_Silica_MH_DHS_f0.7_rv0.1.dat'
# q10 = 'Q_Silica_MH_DHS_f0.7_rv2.0.dat' 
# q11 = 'Q_amorph_c_rv0.1.dat'
# q12 = 'Q_amorph_c_rv1.5.dat'
# cs = ['c_amols_i','c_amoll_i','c_ampys_i','c_ampyl_i','c_fos_i','c_fol_i','c_ens_i','c_enl_i',#'c_amsis_i','c_amsil_i',#'c_amcs_i','c_amcl_i',
    # 'c_amols_o','c_amoll_o','c_ampys_o','c_ampyl_o','c_fos_o','c_fol_o','c_ens_o','c_enl_o']#,'c_amsis_o','c_amsil_o']#,'c_amcs_o','c_amcl_o']
# qs = [q1,q2,q3,q4,q5,q6,q7,q8]*2 #q9,q10,q11,q12]*2

#q1 = 'dustkappa_pure_Fe.dat'
#q2= 'dustkappa_0.05-1um.dat' #'dustkappa_ol-c-mg95_0.1-10um_stitched.dat'  #'dustkappa_pure_Fe.dat' #'dustkappa_0.05-1um.dat' #'dustkappa_default.dat' 
#q3= 'dustkappa_default.dat' #'dustkappa_0.05-1um.dat' #'dustkappa_default.dat' #
#ki1 = 'fe-c_0.10-100.00um.dat' 
#ki2 = 'pyr-mg100_0.20-2.00um.dat'
#ki3 = 'ol-c-mg95_0.77_ens_0.23_-p_0.25_0.20-2.00um.dat' #'ol-c-mg95_0.20-2.00um.dat' #'pyr-mg100_0.20-100.00um.dat' #'ol-mg50_0.20-100.00um.dat' #'pyr-mg70_0.20-500.00um.dat' #'pyr-mg70_0.87_fe-c_0.13_0.05-3000.00um.dat' #'pyr-mg70_0.87_c_0.13_0.05-3000.00um.dat'


'''
km1 = 'fe-c_0.10-100.00um.dat' 
km2 = 'pyr-mg100_0.20-10.00um.dat'
km3 = 'ol-c-mg100_0.20-100.00um.dat'
ko1 = 'fe-c_0.10-100.00um.dat'  #'pyr-mg70_0.87_c_0.13_0.05-3000.00um.dat' #'fe-c_0.10-0.30um.dat'
ko2 = 'pyr-mg100_0.20-10.00um.dat'
ko3 = 'ol-c-mg100_0.20-100.00um.dat'
'''
'''
ki1 = 'dustkappa_pure_Fe.dat'
ki2 = 'dustkappa_default.dat'
ko1 = 'dustkappa_pure_Fe.dat'
ko2 = 'dustkappa_default.dat'
'''
       
#kappas:
# 'pyr-mg70_0.87_fe-c_0.13_0.05-0.30um.dat' not good, too strong sil feat and too small 6-7 um flux
# 'pyr-mg70_0.87_fe-c_0.13_0.05-3.00um.dat' not good, too strong sil feat and too small 6-7 um flux
# in:'fe-c_0.10-0.30um.dat', out: 'pyr-mg70_0.87_fe-c_0.13_0.05-3000.00um.dat', still not good, too large bump at 3-4 um 
# in:'fe-c_0.10-0.30um.dat', out: 'pyr-mg70_0.87_c_0.13_0.05-3000.00um.dat': good fit!

cs =  ['sdens_z%d_%d'%(y+1,x+1) for y in range(N_zones) for x in range(N_dust)]
#cs=['sdens_z1_1','sdens_z2_1','sdens_z3_1']#'sdens_z1_2','sdens_z1_3','sdens_z2_1','sdens_z2_2','sdens_z2_3','sdens_z3_1','sdens_z3_2','sdens_z3_3']#,'c_dustkappa_c']
#cs=['sdens_i1','sdens_i2','sdens_o1','sdens_o2']
#qs=[ki1]*3 #,km1,km2,km3,ko1,ko2,ko3]
qs=ks*N_zones

#available_models: 'TGMdust' 
# pathq ='/Users/jvarga/Data/MATISSE/q-values/'
# 'pathq':'/allegro6/matisse/varga/pro/opacities/'
# 'pathq':'/allegro6/matisse/varga/pro/q-values/'
# 'pathq':'/allegro6/matisse/varga/pro/QVAL/'
# opac_type: 'kappa': kappa [cm^2 g-1], or 'Q': MAC
# 'use_powerlaw_tempgrad': True: use power-law tempgrad, False: calculate the temperture profile self-consistently
# 'parkappa_idx': index of the parametric kappa in the list ks, default: None
model_options = {'nu0':c0/(wl_dic['wavelengths'][4]/1e6),'nu':np.nan,'nx':fitter_settings['img_size_px'],
    'ring_width_ratio':0.067,'fit_SED':fitter_settings['fit_SED'],
    'scale_mas':None, 'calc_profile':False,'pathq':'/allegro6/matisse/varga/pro/q-values/', #binned_MAT_LR-N/full_wlrange/',
    'opac_type':'Q','cs':cs,'qs':qs,'Nzones':N_zones,'use_powerlaw_tempgrad':True,'parkappa_idx':parkappa_idx,
    'V_CF_chi2_weight':3.0,'wl_min':8.0,'wl_max':13.0,'wl_dic':wl_dic,
    #'posterior_range_limits' : [{'param_key':'r_in','min':0.1,'max':math.inf},
    #                            {'param_key':'tau_mp','min':1e-8,'max':1.0},
    #                            {'param_key':'C_i','min':1e-8,'max':0.002}]
    #                            {'param_key':'C_o','min':1e-8,'max':0.0004}
    #                           ]
                            #{'param_key':'r_qjump','min':0.001,'max':2.0},
    #                    [{'param_key':'c_fos_i','min':-math.inf,'max':100.2},
    #                      {'param_key':'c_ens_i','min':-math.inf,'max':2.0},
    #                      {'param_key':'c_enl_i','min':-math.inf,'max':2.0},
    #                      ]
    }

#calculate int(kappa_nu * B_nu(T*) dnu)
fpath = outputdir+'/'+'model_opacities_'+run_name+'.png'
fig, ((ax)) = plt.subplots(1,1,figsize=(10,4))
print(fpath)
model_options['temp_cal_pm'] = []
for i in range(len(qs)):
    input_path = model_options['pathq'] + '/' + qs[i]
    if qs[i] == 'parametric':
        opac_orig = wl_orig*np.nan
    else:
        if model_options['opac_type'] == 'kappa':
            cl=0
            f1 = open(input_path,'r')
            lines = f1.readlines()
            for line in lines:
                if line[0] == '#':
                    cl+=1
            f1.close()
            wl,opac = np.loadtxt(input_path, comments="#", skiprows=cl+2, usecols=(0,1), unpack=True)
        if model_options['opac_type'] == 'Q':
            N,gsize,gdens=np.loadtxt(input_path,max_rows=1)
            wl,Q = np.loadtxt(input_path, comments="#", skiprows=1, usecols=(0,1), unpack=True)
            #convert Q to kappa
            opac = Q* 3.0 / (4.0 * gsize * 1e-4 * gdens)
        
        if i == 0:
            wl_ref = wl
            wl_orig = wl
            opac_orig = opac
        else:
            opac_inter = interpolate.interp1d(wl,opac,kind='linear',fill_value='extrapolate')
            opac_orig = opac_inter(wl_ref)
            wl_orig = wl_ref
        
    if i < N_dust:
        ax.plot(wl_orig,opac_orig,label=qs[i])
            
    #opac_orig = opac_orig*0.0+1.0 #test grey opacity
    nu = c0/(wl_orig*1e-6) # (Hz)
    Bnu_star = B_nu(T_star,nu)
    leftside = (R_star*0.00465047)**2*np.trapz(opac_orig*Bnu_star,x=nu)
    model_options['temp_cal_pm'].append({'leftside':leftside,'wl_orig':wl_orig,
        'nu_orig':nu,'opac_orig':opac_orig,'Bnu_star':Bnu_star})
    #model_options['leftside_in'] = leftside_in
    #model_options['wl_orig'] = wl_orig
    #model_options['nu_orig'] = nu
    #model_options['opac_orig_in'] = opac_orig_in
    #model_options['Bnu_star'] = Bnu_star

ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel(r'Wavelength ($\mu$m)')
ax.set_ylabel(r'$\kappa_\mathrm{abs}$ (cm$^2$ g$^{-1}$)')
ax.set_ylim((1.0,3e4))
ax.set_xlim((0.2,1e3))
plt.legend(loc='best')
plt.tight_layout(pad=0.2)
plt.savefig(fpath)
#plt.savefig(outdir+fname+'.eps',dpi=300)
plt.close(fig)

def calc_cooling_factors(model_options,T_dust,T_star):
    Nc = len(model_options['cs'])
    #print(Nc,T_dust,T_star)
    for j in range(Nc):
        kappa = model_options['temp_cal_pm'][j]['opac_orig']
        nu = model_options['temp_cal_pm'][0]['nu_orig']
        #print(kappa,nu)
        epsilon = np.trapz(kappa*B_nu(T_dust,nu),x=nu) / T_dust**4 / np.trapz(kappa*B_nu(T_star,nu),x=nu) * T_star**4
        print(model_options['qs'][j],": epsilon = %.3f"%epsilon)

print('Dust cooling factors at 1500 K')
calc_cooling_factors(model_options,1500.0,T_star)
print('Dust cooling factors for the dust mix at 1500 K')
T_dust = 1500.0
Nc = len(model_options['cs'])
for kk in range(N_zones):
    kappa = 0.0*model_options['temp_cal_pm'][0]['opac_orig']
    nu = model_options['temp_cal_pm'][0]['nu_orig']
    for j in range(kk*int(Nc/N_zones),(kk+1)*int(Nc/N_zones)):
        kappa += model_options['temp_cal_pm'][j]['opac_orig']* \
            10.0**sdens_dic['sdens_z%d_%d_init'%(kk+1,j%(Nc/N_zones)+1)]
    epsilon = np.trapz(kappa*B_nu(T_dust,nu),x=nu) / T_dust**4 / np.trapz(kappa*B_nu(T_star,nu),x=nu) * T_star**4
    wl_orig = c0/nu/1e-6
    idx = np.logical_and(wl_orig > 1.6, wl_orig < 13.0)
    #print(wl_orig[idx][::2])
    #print(kappa[idx][::2])
    print("Zone %d: epsilon = %.3f"%(kk+1,epsilon))

'''
input_path = model_options['pathq'] + '/' + qs[1]
cl=0
f1 = open(input_path,'r')
lines = f1.readlines()
for line in lines:
    if line[0] == '#':
        cl+=1
f1.close()
wl_orig,opac_orig_out = np.loadtxt(input_path, comments="#", skiprows=cl+2, usecols=(0,1), unpack=True)
leftside_out = (R_star*0.00465047)**2*np.trapz(opac_orig_out*Bnu_star,x=nu)
model_options['leftside_out'] = leftside_out
model_options['opac_orig_out'] = opac_orig_out
'''
'''
input_path = model_options['pathq'] + '/' + qs[1]
cl=0
f1 = open(input_path,'r')
lines = f1.readlines()
for line in lines:
    if line[0] == '#':
        cl+=1
f1.close()
wl_orig_cryst,opac_orig_cryst = np.loadtxt(input_path, comments="#", skiprows=cl+2, usecols=(0,1), unpack=True)
nu = c0/(wl_orig_cryst*1e-6) # (Hz)
Bnu_star = B_nu(T_star,nu)
model_options['wl_cryst'] = wl_orig_cryst
model_options['nu_cryst'] = nu
leftside_cryst = (R_star*0.00465047)**2*np.trapz(opac_orig_cryst*Bnu_star,x=nu)
model_options['leftside_cryst'] = leftside_cryst
model_options['opac_orig_cryst'] = opac_orig_cryst
'''
'''
plt.figure()
plt.plot(wl_orig,opac_orig,marker='.')
plt.plot(wl_orig,Bnu_star)
plt.plot(wl_orig,opac_orig*Bnu_star)
plt.yscale('log')
plt.xscale('log')
plt.show()
'''
#available_models: 'gaussian','ring_asymm, 'ring_symm', 'temp_grad_asymm', 'temp_grad_symm' , 'temp_grad_flared_symm' 
#img_options = {'nu0':c0/(wavelengths[0]/1e6),'nu':np.nan,'nx':img_size_px,'ring_width_ratio':0.1,'scale_mas':None, 'calc_profile':False}

########################### Misc. settings #############################

# limit posterior range, list of dicts