# -*- coding: utf-8 -*-
"""
Created on Thu Jun 23 09:55:18 2022

@author: frede
"""
import numpy as np
from astropy.table import Table

'''

This program calculates the scaling factor of the radiation field G0 relative

to the data of an object or a set of objects radiating for given distances
A file containing my G0 values with the associated distance will then be created

We take the standard Habing value for the interstellar medium 5.6x10^-14 ergs cm-3

between 91.2 and 240nm as normalization factor of G0

The file to be used must contain in column 1 the wavelength in nm,
and in column 2 the intensity per wavelength per sr (in erg cm-2 s-1 nm-1 sr-1).

Each column should not have a name

Example of execution of the program for the star of 10 solar radius HD200775,
for a distance to the star at 20pc :

radiation_field(filename='HD200775_RF.txt', star_radius=10, parsec=20)

If I want to generate a hundred values of G0 :

radiation_field('HD200775_RF.txt', 10, 20, ISRF=False, RF_list=True)

The program will tell me that for 20pc my first value of G0 is not of the order of 1e6,
and will make me enter a new value of distance (in pc) until finding G0 of the order of 1e6

If I want to calculate G0 for the interstellar medium
with an approach described by Habing (1968),

(using Habing file with datas of wavelength and intensity per wavelength per sr) :

radiation_field('habing1968.txt', 1, 1, ISRF=True)

For the interstellar medium, the parameters star_radius, parsec and RF_list

are no longer important;

the part of the program executed will be independent of them

'''

''' Constants '''
h = 6.62607015e-34 #Planck constant in J s
c = 299792458 #Light speed in m s-1

''' Conversions '''
ev = 1.602176634e-19 # 1 ev = 1.602176634e-19 J and value of electron charge

def radiation_field(filename, star_radius , parsec, ISRF=False, RF_list=False):

    """

    ----------

    filename : str,

        name of the .txt containing intensity and wavelength data
    parsec : float,
        distance in parsec from the star
    star_radius : float,
        star radius, in unit of solar radius

    ISRF : bool,

        Interstellar Radiation Field; if true, then the filename is a file for ISRF
        if false, the radiation field of a star is studied
    RF_list: bool,
        If false, the calculation of G0 is done for a given distance,

        if true, it is done for a hundred distances whose input distance
        will correspond to a G0 of the order of 1e6

    returns the value of G0 for a given distance, the distance,
    wavelengths from 91.2nm (13.6ev) to 2 microns, the corresponding intensity

    at a distance d_0, the energy intensity, the energy, the values of ISRF and RF_list

    -------

    """

    ''' association of file data '''

    wave_intensity = Table.read(filename, format='ascii')

    wavelength = wave_intensity['col1'] #in nm
    wavelength_intensity = wave_intensity['col2'] #in erg cm-2 s-1 nm-1 sr-1

    #=> intensity per nm per sr

    ''' useful parameters '''

    radiation_field = []
    distance = []
    energy = ( ( (h/ev) * (c*1e2) ) / (wavelength*1e-7) )[::-1] #in ev
    r = star_radius*7e10 #radius of the star in cm
    d_0 = 3.086e18*parsec #1pc = 3.086e18cm

    ''' energy intensity in erg s-1 cm-2 sr-1 eV-1 '''

    energy_intensity = wavelength_intensity[::-1] * ( (h/ev) * (c*1e9) )/(energy**2)

    ''' indentation '''
    i=0

    if ISRF :

        ''' Radiation Field '''
        intensity_range = (wavelength >= 91.2) & (wavelength <= 240)
        #far UV range, 91.2nm corresponding to 13.6ev
        g_0 = np.trapz( 2*np.pi * wavelength_intensity[intensity_range], wavelength[intensity_range])/(5.6e-14 * c*1e2)
        #g_0 dimensionless
        radiation_field.append(g_0)
        distance.append(0)

    else:

        if RF_list :
            while i < 100:

                ''' distance '''
                if i == 0 :
                    d = d_0
                else:
                    d = 5*i * d_0

                ''' initialization '''
                g_0 = 0

                ''' geometrical dillution '''
                diluted_intensity = wavelength_intensity*(r/d)**2
                #intensity diluted and extincted in erg cm-2 s-1 nm-1 sr-1

                '''energy intensity in erg s-1 cm-2 sr-1 eV-1'''
                energy_intensity = diluted_intensity[::-1] * ( (h/ev) * (c*1e9) )/(energy**2)

                ''' Radiation Field '''
                intensity_range = (wavelength >= 91.2) & (wavelength <= 240)
                g_0 = np.trapz(2 * np.pi * diluted_intensity[intensity_range], wavelength[intensity_range])/(5.6e-14 * c * 1e2)
                distance.append(d)
                radiation_field.append(g_0)

                if (radiation_field[0] >= 2e6 or radiation_field[0] < 1e6):
                    radiation_field = []
                    parsec = float(input('radiation_field[0] = {}, set a value of parsec such that radiation_field[0] is proportional to 1e6 : '.format(g_0),))
                    distance = []
                    d_0 = 3.086e18*parsec
                    i = 0
                else:

                    i = i + 1

        else:

            d = d_0
            ''' geometrical dillution '''
            diluted_intensity = wavelength_intensity*(r/d)**2

            '''energy intensity in erg s-1 cm-2 sr-1 eV-1'''
            energy_intensity = diluted_intensity[::-1] * ( (h/ev) * (c*1e9) )/(energy**2)

            ''' Radiation Field '''
            intensity_range = (wavelength >= 91.2) & (wavelength <= 240)
            g_0 = np.trapz(2 * np.pi * diluted_intensity[intensity_range], wavelength[intensity_range])/(5.6e-14 * c * 1e2)
            distance.append(d)
            radiation_field.append(g_0)

    ''' creation of data files '''

    document_1 = Table([wavelength , wavelength_intensity] , names =('col1','col2'))
    document_1.meta['comments'] = ['wavelength (nm)' , 'intensity (erg s-1 cm-2 nm-1 sr-1)']
    document_1.write('intensity_per_wavelength.txt', format = 'ascii', overwrite=True)

    document_2 = Table([distance , radiation_field], names =('col1', 'col2'))
    document_2.meta['comments'] = ['distance (cm)' , 'radiation field (dimensionless)']
    document_2.write('radiation_field_per_distance.txt', format = 'ascii', overwrite=True)

    return radiation_field[0], distance[0], wavelength, wavelength_intensity, energy_intensity, energy, ISRF, RF_list