#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Oct  7 09:46:46 2020

@author: sacha
"""

import numpy as np
import matplotlib.pyplot as plt
#import matplotlib as mpl
#import coupes_fits as cf
import photoionization_lib as photo_lib
#from astropy.io import fits

# =============================================================================
'''some constants'''
habing=1.6e-3 #habing flux
#Nc=54 #number of carbon atoms of the PAH consedred

#fraction of the remaining energy "beared" by the photo-electron
eV_to_erg=1.6021765e-12 #1 eV in erg 
c=299792458000000000 #light speed in nm
h=4.135667516e-15 #planck constant in eV*sec
h_J=6.626070e-34 #planck constant in J*sec
#
## ============================================================================

'''_________________________________________________________________________'''
'''_____________________________NGC 7023 gas heating model__________________'''
'''_________________________________________________________________________'''
'''physical conditions'''
arcsec=np.array([20,25,30,35,40,45,50,55,60])
d=3.40486e-3*arcsec*3.086e18 #distance to the star in cm for 0.143pc at 42" and d=430pc
T_list=[750,750,750,750,750,750,350,150,150]
Av_list=[0,0,0,0,0,0,0,0.5,1,1.6]#        d=430*arcsec*4.868e-6*3.086e18 #distance to the star in cm
nH_list=np.array([1.5e4,1.5e4,1.5e4,1.5e4,2e4,2e4,2e4,2e4,2e4])
ne_list=1.6e-4*nH_list
#
##carbon locked in PAHs
fc_PAH=np.array([3.9,6,7,7,7,7,7,7,7])*1e-2 #from Berné et al. 2015


wave_range=[91.16,2000] #in nm
#   energy vector

G0_wave_l=[91.16,240] #in nm
G0_wave_E=[h*c/240,h*c/91.16] #in eV
#


'''_________________________________________________________________________'''
'''_____________________________PAH parameters______________________________'''
'''selection of the partition coefficient'''
G=0.46 #usually fixed at 0.5 (Verstraete et al.1990, Tielens (2005))
print('Model started with a partition coefficient of',G)
'''selection of the PAH size'''
Nc_list=np.array([54])
Nc_index=0
'''computation of the photoabsorption crosssection'''
eVN_,eVC_,eVDC_,cross_N_,cross_C_,crossDC_=photo_lib.cross_secs(Nc_list[Nc_index],quiet=True) #/!\ cross sections in Mega Barn => 1e-18 cm2 
'''Estimation of the IPs'''
IP_N=photo_lib.IP_estimate(Nc_list[Nc_index],Z=0)
IP_C=photo_lib.IP_estimate(Nc_list[Nc_index],Z=1)
'''Estimation of the photoionization yield'''

yield_N_=photo_lib.yield_n_to_p(eVN_,IP_N)
yield_C_=photo_lib.yield_p_to_2p(eVC_,IP_C,Nc_list[Nc_index])

'''_________________________________________________________________________'''
'''_____________________________Saving arrays_______________________________'''
nb_point=len(arcsec) #number of modeled locations
gaz_heating_rate_7023=np.zeros([2,nb_point])
#
total_gas_heating_BT=np.zeros(nb_point) # Array of gas heating rate to be compared with Bakes & Tielens 94
total_gas_heating_WD=np.zeros(nb_point) # Array of gas heating rate to be compared with Weingartner & Draine 01

#tableau des taux de chauffage len(Nc_list) la poussiere de chaque etat d'ionisation en chaque points
dust_heating_rate_7023=np.zeros([3,nb_point])

#tableau des taux d'ionisation de chaque etat d'ionisation en chaque points
ionisation_rate_7023=np.zeros([2,nb_point]) #seuls les neutres et cations peuvent etre ionisés

#tableau des taux de recombinaison de chaque etat d'ionisation en chaque points
recombinaison_rate_7023=np.zeros([2,nb_point]) #seuls les cations et dications peuvent se recombiner

#fraction de chaque etat d'ionisation en chaque points
population_fraction_7023=np.zeros([3,nb_point])
ioni_fraction_7023=np.zeros(nb_point)
#efficacité intrinseque de chaque etat d'ionisation en chaque points
intr_efficiency_7023=np.zeros([2,nb_point])
#efficacité intrinseque de chaque etat d'ionisation en chaque points
total_efficiency_7023=np.zeros([nb_point])
#mean photoelectron energy
mean_photelec_energy_7023=np.zeros([2,nb_point])
#total heating rate:
total_gas_heating=np.zeros([nb_point])


#G0 de NGC7023:
G0_list_7023=np.zeros([nb_point])
#gamma array
gamma_list_7023=np.zeros([nb_point])
'''_________________________________________________________________________'''




'''________Application of the model for nb_point points in NGC 7023 NW______'''
for point in range(len(arcsec)):
    T=T_list[point]
    ne=ne_list[point]
    nH=nH_list[point]
    '''====================================================================='''
    '''=======================|radiation field|============================='''
    Fe,E,Fe_l,wave_l=photo_lib.preparation_Fe(d[point],Av_list[point],wave_range=wave_range,quiet=True) #Fe en erg/sec/cm2/eV/sr
#        Fe=Fe_*2*np.pi
    G0=np.trapz(2*np.pi*Fe[(E>G0_wave_E[0]) & (E<G0_wave_E[1])],E[(E>G0_wave_E[0]) & (E<G0_wave_E[1])])/1.6e-3
    G0_list_7023[point]=G0.copy()
    gamma_list_7023[point]=np.sqrt(T)*G0/ne
    print('G0={} and gamma={} at {}"'.format(G0,gamma_list_7023[point],arcsec[point]))

    '''====================================================================='''
    '''=========================|ionization cross-section building|========='''
    '''====================================================================='''
    ''' adaptating the cross section from 0 to 13.6eV''' 
    crossN=np.interp(E,eVN_,cross_N_)*1e-18 #in cm2/Carbon (from Mb/C to cm2/C)
    crossC=np.interp(E,eVC_,cross_C_)*1e-18 #in cm2/Carbon (from Mb/C to cm2/C)
    crossDC=np.interp(E,eVDC_,crossDC_)*1e-18 #in cm2/Carbon (from Mb/C to cm2/C)
    yield_N=np.interp(E,eVN_,yield_N_) 
    yield_C=np.interp(E,eVC_,yield_C_) 
    
    
    ioni_cross_N=crossN.copy()*yield_N.copy()     #in cm2/C (from Mb/C to cm2/C)
    ioni_cross_C=crossC.copy()*yield_C.copy()     #in cm2/C (from Mb/C to cm2/C)
    
    
    '''====================================================================='''
    '''=====================dust and gas heating calculation================'''
    '''====================================================================='''
    E_range_N=np.where(E>IP_N)[0]
    E_range_C=np.where(E>IP_C)[0]
#        E_range_DC=np.where(E>IP_DC)[0]
    
    '''photo absorption of the neutrals, cations and dications '''
    photo_absorption_N=Fe*crossN*2*np.pi #erg/sec/eV /!\ the 2*np.pi is the solid angle considered =>  the RF comes from the star only
    photo_absorption_C=Fe*crossC*2*np.pi #erg/sec/eV
    photo_absorption_DC=Fe*crossDC*2*np.pi #erg/sec/eV
        
    '''heating rate of the molecule itself'''
    chauffe_dust_N=np.trapz(photo_absorption_N,E) #erg/sec
    chauffe_dust_C=np.trapz(photo_absorption_C,E) #erg/sec
    chauffe_dust_DC=np.trapz(photo_absorption_DC,E) #erg/sec
 
    
    #saving of the dust heating rate
    dust_heating_rate_7023[0,point]=chauffe_dust_N.copy()
    dust_heating_rate_7023[1,point]=chauffe_dust_C.copy()
    dust_heating_rate_7023[2,point]=chauffe_dust_DC.copy()
    
    
    ''' power density absorbed for ionization '''
    ionisation_absorption_N=yield_N*photo_absorption_N #erg/sec/eV
    ionisation_absorption_C=yield_C*photo_absorption_C #erg/sec/eV
#    ionisation_absorption_DC=yield_2p_to_3p*photo_absorption_DC #erg/sec/eV
    
    '''number of ionizations'''
    number_ionisation_absorption_N=ionisation_absorption_N/(E*eV_to_erg) # ioni/sec/eV
    number_ionisation_absorption_C=ionisation_absorption_C/(E*eV_to_erg) # ioni/sec/eV
#    number_ionisation_absorption_DC=ionisation_absorption_DC/(E*eV_to_erg) # ioni/sec/eV
    
    '''photoionization rate'''
    kpe_N=np.trapz(number_ionisation_absorption_N[E_range_N],(E-IP_N)[E_range_N])
    kpe_C=np.trapz(number_ionisation_absorption_C[E_range_C],(E-IP_C)[E_range_C])
#        kpe_DC=np.trapz(number_ionisation_absorption_DC[E_range_DC],(E-IP_DC)[E_range_DC])

    '''mean gas heating energy'''
    mean_photelec_energy_7023[0,point]=1/kpe_N*np.trapz(G*(E[E_range_N]-IP_N)*number_ionisation_absorption_N[E_range_N],(E[E_range_N]-IP_N)) #from Verstraete et al. 90
    mean_photelec_energy_7023[1,point]=1/kpe_C*np.trapz(G*(E[E_range_C]-IP_C)*number_ionisation_absorption_C[E_range_C],(E[E_range_C]-IP_C))
    
    #saving of the ionization rate
    ionisation_rate_7023[0,point]=kpe_N.copy()
    ionisation_rate_7023[1,point]=kpe_C.copy()
#        ionisation_rate[2,mol,point]=kpe_DC.copy()

    '''spectrum of the gas heating per charge state'''
    chauffage_du_gaz_par_eV_N=G*(E-IP_N)*number_ionisation_absorption_N*eV_to_erg #erg/sec/eV
    chauffage_du_gaz_par_eV_C=G*(E-IP_C)*number_ionisation_absorption_C*eV_to_erg #erg/sec/eV
#        chauffage_du_gaz_par_eV_DC=G*(E-IP_DC)*number_ionisation_absorption_DC*eV_to_erg #erg/sec/eV
    
    '''gaz heating rate per charge state'''
    chauffe_gas_neutre=np.trapz(chauffage_du_gaz_par_eV_N[E_range_N],(E[E_range_N]-IP_N)) #erg/sec/molecule
    chauffe_gas_cation=np.trapz(chauffage_du_gaz_par_eV_C[E_range_C],(E[E_range_C]-IP_C)) #erg/sec/molecule
#        chauffe_gas_dication=np.trapz(chauffage_du_gaz_par_eV_DC[E_range_DC],G*(E[E_range_DC]-IP_DC))
    
    #saving of the gas heating rate
    gaz_heating_rate_7023[0,point]=chauffe_gas_neutre.copy()
    gaz_heating_rate_7023[1,point]=chauffe_gas_neutre.copy()

      
    '''====================================================================='''    
    '''======================= population fractions ========================'''
    '''======================= in a 3-levels model ========================='''
    '''====================================================================='''
    
    '''======================= recombinaison ========================'''
    krec_neutral=ne*1.28e-10*Nc_list[Nc_index]*np.sqrt(T)*(1+(1.85e5)/(np.sqrt(Nc_list[Nc_index])*T))
    krec_cation=ne*1.28e-10*Nc_list[Nc_index]*np.sqrt(T)*(1+(2*1.85e5)/(np.sqrt(Nc_list[Nc_index])*T))
#     krec_dication=ne*1.28e-10*Nc_list[Nc_index]*np.sqrt(T)*(1+(3*1.85e5)/(np.sqrt(Nc_list[Nc_index])*T))

    #saving of the recombinaison rate
    recombinaison_rate_7023[0,point]=krec_neutral.copy()
    recombinaison_rate_7023[1,point]=krec_cation.copy()
#                recombinaison_rate[temp,radfield,2,Nc_index,point]=krec_dication.copy()
    
    '''population fraction computation'''
    #cations
    fp3l=1/(1+krec_neutral/kpe_N+kpe_C/krec_cation)
    #neutrals
    f03l=(1-fp3l*kpe_C/krec_cation)/(1+kpe_N/krec_neutral)
    #dications
    fpp3l=(1-f03l)/(1+krec_cation/kpe_C)

    #saving of the fraction of each population
    population_fraction_7023[0,point]=f03l.copy()
    population_fraction_7023[1,point]=fp3l.copy()
    population_fraction_7023[2,point]=fpp3l.copy()
    '''ionization fraction'''
    ioni_fraction_7023[point]=fpp3l.copy()+fp3l.copy()

    '''======================= heating efficiencies ========================'''
    '''intrinsic'''
    intr_efficiency_N=chauffe_gas_neutre/chauffe_dust_N
    intr_efficiency_C=chauffe_gas_cation/chauffe_dust_C
    
    #saving of the intrinsic heating efficiency of each population
    intr_efficiency_7023[0,point]=intr_efficiency_N.copy()
    intr_efficiency_7023[1,point]=intr_efficiency_C.copy()
    
    '''mixed 3-levels model'''
    total_gas_heating[point]=(f03l*chauffe_gas_neutre+fp3l*chauffe_gas_cation) #erg/sec/(molecule de taille Nc)
    
    mixed_efficiency_3l=(f03l*chauffe_gas_neutre+fp3l*chauffe_gas_cation)/(f03l*chauffe_dust_N+\
                        fp3l*chauffe_dust_C+fpp3l*chauffe_dust_DC)
    total_efficiency_7023[point]=mixed_efficiency_3l.copy()

    '''gas heating'''
    #saving of the total heating efficiency
    #
    total_gas_heating_WD[point]=(f03l*chauffe_gas_neutre+fp3l*chauffe_gas_cation)/(G0*nH)
    total_gas_heating_BT[point]=(f03l*chauffe_gas_neutre+fp3l*chauffe_gas_cation)*fc_PAH[point]/Nc_list[Nc_index]*1.6e-4

'''display of the principal parameters'''
''' pop_fraction with the distance ''' 
fig0,ax1=plt.subplots(figsize=(5.90,3.93))
ax1.set_xlabel('Distance fron the star (")')
ax1.set_ylabel('population fraction')
ax1.plot(arcsec,population_fraction_7023[0,:],color = 'green',label='Z=0')
ax1.plot(arcsec,population_fraction_7023[1,:],color = 'blue',label='Z=1')
ax1.plot(arcsec,population_fraction_7023[2,:],color = 'red',label='Z=2')
plt.legend()            
                
'''comparaison with observation'''
if True:
    '''photoionization yield''' 
    IP_Npet=photo_lib.IP_estimate(32,0)
    IP_Nmoy=photo_lib.IP_estimate(44,0)
    IP_Ngrand=photo_lib.IP_estimate(64,0)
    
    IP_Cpet=photo_lib.IP_estimate(32,1)
    IP_Cmoy=photo_lib.IP_estimate(44,1)
    IP_Cgrand=photo_lib.IP_estimate(64,1)
    
    yield_petitN=photo_lib.yield_n_to_p(eVN_,IP_Npet)
    yield_moyN=photo_lib.yield_n_to_p(eVN_,IP_Nmoy)
    yield_grandN=photo_lib.yield_n_to_p(eVN_,IP_Ngrand)
    
    yield_petitC=photo_lib.yield_p_to_2p(eVC_,IP_Cpet,32)
    yield_moyC=photo_lib.yield_p_to_2p(eVC_,IP_Cmoy,44)
    yield_grandC=photo_lib.yield_p_to_2p(eVC_,IP_Cgrand,66)

    plt.figure()
    plt.plot(eVN_,yield_petitN,'darkblue',label='Z=0->1')
    plt.plot(eVN_,yield_moyN,'cornflowerblue')
    plt.plot(eVN_,yield_grandN,'darkred')
    plt.text(4.1,0.70,'N$_C$=32',color='darkblue')
    plt.text(4.1,0.65,'N$_C$=44',color='cornflowerblue')
    plt.text(4.1,0.60,'N$_C$=66',color='darkred')
    plt.plot(eVN_,yield_petitC,'darkblue',ls='--',label='Z=1->2')
    plt.plot(eVN_,yield_moyC,'cornflowerblue',ls='--')
    plt.plot(eVN_,yield_grandC,'darkred',ls='--')
    plt.xlabel('Énergie de photon (eV)')
    plt.ylabel('Y(E,Z)')
    plt.xlim([4,13.6])
    plt.legend()
    '''====================================================================='''    

    
    '''ioni and photo cross sections'''
    eVN_,eVC_,eVDC_,cross_N_petit,cross_C_petit,crossDC_petit=photo_lib.cross_secs(32,quiet=True) #/!\ cross sections en Mega Barn => 1e-18 cm2 
    cross_N_petit=cross_N_petit/32*1e-18
    cross_C_petit=cross_C_petit/32*1e-18
    crossDC_petit=crossDC_petit/32*1e-18
    ionicross_N_petit=cross_N_petit*yield_petitN
    ionicross_C_petit=cross_C_petit*yield_petitC

    
    eVN_,eVC_,eVDC_,cross_N_moy,cross_C_moy,crossDC_moy=photo_lib.cross_secs(44,quiet=True) #/!\ cross sections en Mega Barn => 1e-18 cm2 
    cross_N_moy=cross_N_moy/44*1e-18
    cross_C_moy=cross_C_moy/44*1e-18
    crossDC_moy=crossDC_moy/44*1e-18
    ionicross_N_moy=cross_N_moy*yield_moyN
    ionicross_C_moy=cross_C_moy*yield_moyC
    
    eVN_,eVC_,eVDC_,cross_N_grand,cross_C_grand,crossDC_grand=photo_lib.cross_secs(66,quiet=True) #/!\ cross sections en Mega Barn => 1e-18 cm2 
    cross_N_grand=cross_N_grand/66*1e-18
    cross_C_grand=cross_C_grand/66*1e-18
    crossDC_grand=crossDC_grand/66*1e-18
    ionicross_N_grand=cross_N_grand*yield_grandN
    ionicross_C_grand=cross_C_grand*yield_grandC
    
    fig=plt.figure(figsize=(3.97,4.17))
#    fig=plt.figure()
    ax1 = fig.add_subplot(311)
    ax1.plot(eVN_,cross_N_petit*1e17,'darkblue',label='$\sigma(E,Z)$')
    ax1.plot(eVN_,cross_N_moy*1e17,'cornflowerblue')
    ax1.plot(eVN_,cross_N_grand*1e17,'darkred')    
    ax1.plot(eVN_,ionicross_N_petit*1e17,'darkblue',ls='--',label='$\sigma_{ion}(E)/N_C$')
    ax1.plot(eVC_,ionicross_N_moy*1e17,'cornflowerblue',ls='--')
    ax1.plot(eVDC_,ionicross_N_grand*1e17,'darkred',ls='--')    
#    plt.legend()

    
    ax2 = fig.add_subplot(312)
    ax2.plot(eVN_,cross_C_petit*1e17,'darkblue')
    ax2.plot(eVN_,cross_C_moy*1e17,'cornflowerblue')
    ax2.plot(eVN_,cross_C_grand*1e17,'darkred')    
    ax2.plot(eVN_,ionicross_C_petit*1e17,'darkblue',ls='--')
    ax2.plot(eVC_,ionicross_C_moy*1e17,'cornflowerblue',ls='--')
    ax2.plot(eVDC_,ionicross_C_grand*1e17,'darkred',ls='--')    
    
    ax3 = fig.add_subplot(313)
    ax3.plot(eVN_,crossDC_petit*1e17,'darkblue')
    ax3.plot(eVN_,crossDC_moy*1e17,'cornflowerblue')
    ax3.plot(eVN_,crossDC_grand*1e17,'darkred')    

    ax1.tick_params(labelbottom='off',direction='inout')
    ax2.tick_params(labelbottom='off',direction='inout')
    ax3.tick_params(direction='inout')
    ax1.xaxis.set_ticks_position('both')
    ax2.xaxis.set_ticks_position('both')
    ax3.xaxis.set_ticks_position('both')

#    ax1.set_xlabel('Énergie de photon (eV)')
    ax1.text(8,1,'Z=0')
    ax2.text(8,1,'Z=1')
    ax3.text(8,1,'Z=2')
    ax1.text(0,1.2,'- $\sigma_{abs}$')
    ax1.text(0,0.9,'-- $\sigma_{ion}$')
    
#    ax2.set_xlabel('Énergie de photon (eV)')
    ax3.set_xlabel('Photon energy (eV)')
    ax2.set_ylabel('Cross sections (10$^{-17}$cm$^2$/C)')
    ax2.text(0,1.15,'N$_C$=32',color='darkblue')
    ax2.text(0,0.9,'N$_C$=44',color='cornflowerblue')
    ax2.text(0,0.65,'N$_C$=66',color='darkred')
    fig.subplots_adjust(top=0.96,bottom=0.11,left=0.135,right=0.985,hspace=0.0,wspace=0.2)
    '''====================================================================='''    

    '''====================================================================='''    
    '''=====================cut in the observations========================='''    
    import coupes_fits as cf #personnal librairy used for the cuts

    #RA/DEC coordinates in deg of the extremes points of the 3 cuts 
    #starting point 
    coord1_1=[315.3945810,68.1785188]
    coord1_2=[315.3797831,68.1748496] 
    coord1_3=[315.3644449,68.1715722]
    #ending point
    coord2_2=[315.39375,68.1667]
    coord2_3=[315.3909639,68.1662128]
    coord2_1=[315.3978010,68.1687543]
    #HD200775 coord.
    coordstar=[315.4038,68.1632]
    #list of coordinates to be used by the cut function
    coords=[[coord1_1,coord1_2,coord1_3],[coord2_1,coord2_2,coord2_3],coordstar]

    #files to consider
    #PE efficiency, neutral and cation integrated intensity
    filename_PE='./data/PE_PAH_h2_PAH_map_publication.fits'
    filename_PAH0='./data/N7023_neutrals_map_intLR_11_reproj.fits'
    filename_PAHp='./data/N7023_cations_map_intLR_11_reproj.fits'    
    #gas emission
    filename_CII='./data/CII157line_BS.fits'
    filename_OI63='./data/OI63line_BS.fits'
    filename_OI145='./data/OI145line_BS.fits'
    filename_AIB='./data/AIB_int_map_publication.fits'
    filename_cont='./data/N7023_int_cont_map.fits'
    
    
    model_output=np.zeros([2,len(arcsec)])
    model_output[0,:]=arcsec.copy()
    model_output[1,:]=total_efficiency_7023.copy()
    
    #cuts, see the documentation and comments of the cf.cut_fits
    #Pe efficiency and neutral and cation integrated intensity
    all_data_average_PE,all_error_PE,all_distance_average_PE=cf.cut_fits(filename_PE,coords,plot_all_cut=True,average=True,model_output=model_output,fig_index='PE efficiency')
    all_data_averagePAH0,all_error_PAH0,all_distance_average_PAH0=cf.cut_fits(filename_PAH0,coords,plot_all_cut=True,average=True,model_output=[],fig_index='PAH0')
    all_data_averagePAHp,all_error_PAHp,all_distance_average_PAHp=cf.cut_fits(filename_PAHp,coords,plot_all_cut=True,average=True,model_output=[],fig_index='PAH+')
    #gas emission
    all_data_average_CII,all_error_CII,all_distance_average_CII=cf.cut_fits(filename_CII,coords,plot_all_cut=True,average=True,model_output=[],fig_index='C+')
    all_data_average_OI63,all_error_OI63,all_distance_average_OI63=cf.cut_fits(filename_OI63,coords,plot_all_cut=True,average=True,model_output=[],fig_index='OI63')
    all_data_average_OI145,all_error_OI145,all_distance_average_OI145=cf.cut_fits(filename_OI145,coords,plot_all_cut=True,average=True,model_output=[],fig_index='OI145')
    all_data_average_AIB,all_error_AIB,all_distance_average_AIB=cf.cut_fits(filename_AIB,coords,plot_all_cut=True,average=True,model_output=[],fig_index='AIB')
    all_data_average_cont,all_error_cont,all_distance_average_cont=cf.cut_fits(filename_cont,coords,plot_all_cut=True,average=True,model_output=[],fig_index='cont')


#
    #
    
    plt.figure(figsize=(5.90,3.93))
    '''verify that you used Nc=54 to be comparable to the following computed PE efficiencies'''
    #PE efficiencies computed for Nc=54 and several partition coefficient
    pe_eff_G04=np.array([0.00613149, 0.00845175, 0.01065419, 0.0126501 , 0.01656323, 0.0176872 , 0.02070102, 0.02283699, 0.02091324])
    pe_eff_G035=np.array([0.00536505, 0.00739529, 0.00932241, 0.01106884, 0.01449282,  0.0154763 , 0.0181134 , 0.01998237, 0.01829909])    
    pe_eff_G05=np.array([0.00766436, 0.01056469, 0.01331773, 0.01581263, 0.02070403, 0.022109  , 0.02587628, 0.02854624, 0.02614155])
    pe_eff_G055=np.array([0.00843079, 0.01162116, 0.01464951, 0.01739389, 0.02277444, 0.0243199 , 0.02846391, 0.03140086, 0.02875571])

    #plot of the computed PAH PE efficiency
    plt.semilogy(population_fraction_7023[1,:]/(population_fraction_7023[0,:]+population_fraction_7023[1,:]),total_efficiency_7023,'x',color='red',label='modèle, G={}'.format(G))
    
    #observations: cuts of the observed ``reduced'' PE efficiency in function of the cation fraction
    plt.errorbar(all_data_averagePAHp[0]/(all_data_averagePAHp[0]+all_data_averagePAH0[0]),all_data_average_PE[0],yerr=all_error_PE[0],label='coupe {}'.format(0))
    plt.errorbar(all_data_averagePAHp[1]/(all_data_averagePAHp[1]+all_data_averagePAH0[1]),all_data_average_PE[1],yerr=all_error_PE[1],label='coupe {}'.format(1))
    plt.errorbar(all_data_averagePAHp[2]/(all_data_averagePAHp[2]+all_data_averagePAH0[2]),all_data_average_PE[2],yerr=all_error_PE[2],label='coupe {}'.format(2))
    #for comparaison to the computed pe_eff_GXX
    plt.plot(population_fraction_7023[1,:]/(population_fraction_7023[0,:]+population_fraction_7023[1,:]),pe_eff_G035,'v',label='model, G=0.35')
    plt.plot(population_fraction_7023[1,:]/(population_fraction_7023[0,:]+population_fraction_7023[1,:]),pe_eff_G04,'1',label='model, G=0.4')
#    plt.plot(population_fraction_7023[1,:]/(population_fraction_7023[0,:]+population_fraction_7023[1,:]),pe_eff_G05,'d',label='modèle, G=0.5')
    plt.plot(population_fraction_7023[1,:]/(population_fraction_7023[0,:]+population_fraction_7023[1,:]),pe_eff_G055,'*',label='model, G=0.55')
    plt.legend()
#    plt.axvline(0.5,lw=0.8,ls='--',color='grey')
    plt.xlabel('Cation fraction')        
    plt.ylabel('$\epsilon$')
    plt.title('N$_C$={}'.format(54))
    


    plt.figure('gas emission')
    plt.semilogy(all_distance_average_PAHp[0],all_data_average_CII[0],label='[CII] 158$\mu$m')
    plt.plot(all_distance_average_PAHp[0],all_data_average_OI63[0],label='[0I] 63$\mu$m')
    plt.plot(all_distance_average_PAHp[0],all_data_average_OI145[0],label='[OI] 145$\mu$m')
    plt.plot(all_distance_average_PAHp[0],all_data_average_AIB[0],label='AIB')
    plt.xlabel('Distance (")')
    plt.ylabel('Brightness (W/m2/sr)')
    plt.legend()

    filename_PAHp_IRSres='./data/N7023_ionization_fraction.fits' #ionization fraction map at the IRS spatial resolution
    all_data_averagePAHp_IRSres,all_error_PAHp_IRSres,all_distance_average_PAHp_IRSres=cf.cut_fits(filename_PAHp_IRSres,coords,plot_all_cut=True,average=True,model_output=[],fig_index='PAH+ IRS')

    #computed ionization fraction for different sizes
    ioni_fraction_NC32=np.array([0.79313524, 0.69150685, 0.59621668, 0.51204338, 0.36650701, 0.31117851, 0.19870353, 0.04037594, 0.01136714])
    ioni_fraction_NC44=np.array([0.85018817, 0.7634044,  0.67614262, 0.59450648, 0.4432105,  0.38221601, 0.25239059, 0.05409562, 0.01544414])
    ioni_fraction_NC66=np.array([0.8855932,  0.81159114, 0.73315735, 0.6564234,  0.50591308, 0.4421752,  0.30113321, 0.06761376, 0.01948429])
    #computed PE on PAHs heating efficiency for different sizes
    PE_eff_NC32_G05=np.array([0.00591002, 0.01120707, 0.01392701, 0.01626842, 0.02021422, 0.02168575, 0.02463521, 0.02668083, 0.02525464])
    PE_eff_NC44_G05=np.array([0.00589942, 0.01120773, 0.0140879 , 0.01667006, 0.0212639 , 0.02306228, 0.02680859, 0.0297366 , 0.02875639])
    PE_eff_NC66_G05=np.array([0.00542426, 0.0103825 , 0.01318074, 0.01577018, 0.02059159, 0.02256117, 0.02680796, 0.03037893, 0.02954722])
    #computed PE on PAHs heating efficiency for Nc=66 and different partition coefficients
    PE_eff_NC66_G06=np.array([0.00650911, 0.012459,   0.01581688, 0.01892422, 0.02470991, 0.02707341, 0.03216955, 0.03645472, 0.03545667])
    PE_eff_NC66_G04=np.array([0.00433941, 0.008306,   0.01054459, 0.01261615, 0.01647328, 0.01804894, 0.02144637, 0.02430315, 0.02363778])
    
    
#    print('$\epsilon_{PAH}$',total_efficiency_7023)
    print('ioni fraction, Nc={}, G={}'.format(Nc_list[0],G))
    
    #comparaison computed ionization fraction VS observed ionization fraction
    plt.figure(figsize=(5.90,3.93))
    ls=['-','--','-.']
    plt.xlabel('Distance (")')
    plt.ylabel("$f\,_{ion}^{Obs}$,$f\,_{ion}^{Mod}$")
    plt.errorbar(all_distance_average_PAHp_IRSres[0],all_data_averagePAHp_IRSres[0],yerr=all_error_PAHp_IRSres[0],label='coupe 1')
    plt.errorbar(all_distance_average_PAHp_IRSres[1],all_data_averagePAHp_IRSres[1],yerr=all_error_PAHp_IRSres[1],label='coupe 2')
    plt.errorbar(all_distance_average_PAHp_IRSres[2],all_data_averagePAHp_IRSres[2],yerr=all_error_PAHp_IRSres[2],label='coupe 3')
#    plt.plot(arcsec,population_fraction_7023[1,:]/(population_fraction_7023[0,:]+population_fraction_7023[1,:]),'x',color='red',label='modèle PAH$^+$')
    plt.plot(arcsec,ioni_fraction_NC66,'x',color='darkred',label='N$_C$=66')
    plt.plot(arcsec,ioni_fraction_NC44,'x',color='cornflowerblue',label='N$_C$=44')
    plt.plot(arcsec,ioni_fraction_NC32,'x',color='darkblue',label='N$_C$=32')
    plt.subplots_adjust(top=0.984,bottom=0.128,left=0.093,right=0.992,hspace=0.2,wspace=0.2)    
    plt.axhline(0.5,ls='--',alpha=0.5)
    plt.legend(loc=3)
    
    
###############
    #plot of the heating efficiencies from WD01 and BT94
    #from WD01
    size_WD01=[4,5,6,7,8,9,10] #in Angstroms
    gamma_WD01=[1e3,1e4,1e5]
    PE_eff_WD01_gam1e3=np.array([0.05,0.06,0.065,0.065,0.065,0.065,0.065])
    PE_eff_WD01_gam1e4=np.array([0.019,0.019,0.019,0.018,0.017,0.016,0.015])
    PE_eff_WD01_gam1e5=np.array([0.0023,0.0021,0.0019,0.0018,0.0016,0.0015,0.0014])
    PE_eff_WD01=np.zeros([7,3])
    PE_eff_WD01[:,0]=PE_eff_WD01_gam1e3.copy()
    PE_eff_WD01[:,1]=PE_eff_WD01_gam1e4.copy()
    PE_eff_WD01[:,2]=PE_eff_WD01_gam1e5.copy()
    
    #from BT94
    #/!\ the size variation is not considered in BT94 formula
    epsilonBT94=(4.87e-2)/(1+(4e-3)*(gamma_list_7023)**0.73)+(3.65e-2*(100/1e4)**0.7)/(1+2e-4*(gamma_list_7023)) 

    plt.figure()
    #plot of WD01
    plt.loglog(gamma_WD01,PE_eff_WD01[0,:],ls='--',color='cornflowerblue',label='WD01, 4$\AA$')
#    plt.plot(gamma_WD01,PE_eff_WD01[1,:],'*')
#    plt.plot(gamma_WD01,PE_eff_WD01[2,:],'*')
    plt.plot(gamma_WD01,PE_eff_WD01[3,:],ls='--',color='darkred',label='WD01, 7$\AA$')
#    plt.plot(gamma_WD01,PE_eff_WD01[4,:],'*')
#    plt.plot(gamma_WD01,PE_eff_WD01[5,:],'*')
    plt.plot(gamma_WD01,PE_eff_WD01[6,:],ls='--',color='darkblue',label='WD01, 10$\AA$')
    #plot of BT94    
    plt.plot(gamma_list_7023,epsilonBT94,ls=':',color='black',label='BT94')
       
    #plot of this code
    plt.plot(gamma_list_7023,PE_eff_NC66_G05,color='darkred',label='modele PAH, Nc=66')
    plt.plot(gamma_list_7023,PE_eff_NC44_G05,color='cornflowerblue',label='modele PAH, Nc=44')
    plt.plot(gamma_list_7023,PE_eff_NC32_G05,color='red',label='modele PAH, Nc=32')
    
    plt.xlabel('$\gamma$ (K$^{1/2}$cm$^3$)')
    plt.ylabel('$\epsilon$')
    plt.legend()
#######    
    
    
    '''introduction of the alpha parameter, i.e. the contribution of PAHs in the gas heating'''
    
    
    alpha=1 #if 1, you consider that only PAH contribute to the gas heating. Good results with alpha=0.8
    plt.figure(figsize=(5.90,3.93))
    plt.errorbar(all_data_averagePAHp[0]/(all_data_averagePAHp[0]+all_data_averagePAH0[0]),all_data_average_PE[0]*alpha,yerr=all_error_PE[0],label='coupe 1'.format(0))
    plt.errorbar(all_data_averagePAHp[1]/(all_data_averagePAHp[1]+all_data_averagePAH0[1]),all_data_average_PE[1]*alpha,yerr=all_error_PE[1],label='coupe 2'.format(1))
    plt.errorbar(all_data_averagePAHp[2]/(all_data_averagePAHp[2]+all_data_averagePAH0[2]),all_data_average_PE[2]*alpha,yerr=all_error_PE[2],label='coupe 3'.format(2))

#    plt.plot(population_fraction_7023[1,:]/(population_fraction_7023[0,:]+population_fraction_7023[1,:]),total_efficiency_7023,'x',color='red',label='$\epsilon_{th}$')
    plt.plot(ioni_fraction_NC66,PE_eff_NC66_G05,'x',color='red',label='Grands')
#    plt.plot(ioni_fraction_NC66,PE_eff_NC66_G05,'x',color='orange',label='G=0.5')
#    plt.plot(ioni_fraction_NC66,PE_eff_NC66_G04,'x',color='blue',label='G=0.4')
    plt.plot(ioni_fraction_NC44,PE_eff_NC44_G05,'x',color='orange',label='Moyens')
    plt.plot(ioni_fraction_NC32,PE_eff_NC32_G05,'x',color='blue',label='Petits')
    plt.subplots_adjust(top=0.984,bottom=0.128,left=0.108,right=0.992,hspace=0.2,wspace=0.2)    
    plt.legend()
    plt.ylabel('$\epsilon\,_{PAH}$, $\epsilon\,_{r}$')
    plt.xlabel("$f\,_{ion}^{Mod}$,$f\,_{ion}^{Obs}$")
    
  
    
    plt.figure()
    plt.loglog(gamma_list_7023,total_gas_heating_WD,label='Nc=66, 8$\AA$')
    plt.ylabel('$\Gamma/G_0n_H$')
    plt.xlabel('$\gamma$ (T$^{1/2}$cm$^3$)')
    
    plt.figure()
    plt.loglog(gamma_list_7023,total_gas_heating_BT,label='Nc=66, 8$\AA$')
    plt.ylabel('taux de chauffage (erg/sec/H)')
    plt.xlabel('$\gamma$ (T$^{1/2}$cm$^3$)')

if False:  
    #T=750
    #G0=2600
    n=7e3
    ne=1.6e-4*n
    nb_temp=50
    T=np.logspace(1,4,nb_temp)
    G0=np.logspace(2,4,nb_temp)
    for i in range(len(G0_list_7023)):
        epsilon=(4.87e-2)/(1+(4e-3)*(G0*np.sqrt(T)/ne)**0.73)+(3.65e-2*(100/1e4)**0.7)/(1+2e-4*(G0_list_7023[i]*np.sqrt(T)/ne))
        chauffageBT=1e-24*epsilon*nH_list[0]*G0
        #chauffage=total_gas_heating*0.1*1.6e-4/Nc_list[0]*n
        refroidissement=2.5e-29*n**2*T**(2/3)*np.exp(-228/T)+8e-27*n**2*np.exp(-92/T)
        plt.figure()
        plt.loglog(T,chauffageBT,'r',label='heating rate')
        plt.plot(T,refroidissement,'b',label='cooling rate')
        plt.legend()

#    fig=plt.figure()
  
    
    nb_temp=1000
    T=np.logspace(1,4,nb_temp)
    refroidissement=np.zeros([len(arcsec),nb_temp])
    
    for j in range(len(arcsec)):
        refroidissement[j,:]=8e-27*nH_list[j]**2*np.exp(-92/T)
        chauffage=total_gas_heating[j]*fc_PAH[point]*2.7e-4/Nc_list[Nc_index]*nH_list[j]
        plt.figure('distance={}'.format(arcsec[j]))
        plt.loglog(T,refroidissement[j,:],'r',label='cooling rate')
        plt.axhline(chauffage,ls='--',label='heating rate')
        plt.legend()
        
heat_rate_7023=total_gas_heating*fc_PAH[point]*2.7e-4/Nc_list[Nc_index]*1e-7 # heating in watts / H