import numpy as np
import matplotlib.pyplot as plt

n_sigma = 30  #number of sigma-values to try
n_flux = 30  #number of flux-min values to try
k=1.38064852*10**(-23)
temp_min = 150.  #minimum temp. of gases
temp_max = 450.  #max temp. of gases


def F_model(wavelength,sigma,F_min,lines):
    F_model= #fill in model here
    return F_model

def Mm(lines): #masses of molecules in kg
    if 630==lines or 690==lines or 760==lines:
        return *1.67e-27
    if 720==lines or 820==lines or 940==lines:
        return *1.67e-27
    if 1400==lines or 1600==lines:
        return *1.67e-27
    if 1660==lines or 2200==lines:
        return *1.67e-27
    if 2340==lines:
        return *1.67e-27
    if 2870==lines:
        return *1.67e-27

def Mn(lines): #names of molecules
    if 630==lines or 690==lines or 760==lines:
        return "Oxygen"
    if 720==lines or 820==lines or 940==lines:
        return "Water vapour"
    if 1400==lines or 1600==lines:
        return "Carbon dioxide"
    if 1660==lines or 2200==lines:
        return "Methane"
    if 2340==lines:
        return "Carbon monoxide"
    if 2870==lines:
        return "Nitrous oxide"

def absorption_line(k=k):
    absorption=np.loadtxt("")  #load .txt-spectrum file
    noise=np.loadtxt("sigma_noise.txt")
    absorption_lines=[630,690,760,720,820,940,1400,1600,1660,2200,2340,2870]
    
    for lines in absorption_lines:
        print('line: ',lines, Mn(lines))
        dopplermax =   #max doppler shift
        fwhm_min = #smallest possible fwhm
        fwhm_max =  #largest possible fwhm
        sigma_min = fwhm_min/ #smallest possible sigma
        sigma_max = fwhm_max/ #largest possible sigma
        sigma = np.linspace(sigma_min, sigma_max ,n_sigma) #array of sigma values to try
        lambda_min = lines - dopplermax - sigma_max*find_a_reasonable_number  #first lambda to try for given line, make sure to capture the whole line
        lambda_max = lines + dopplermax + sigma_max*find_a_reasonable_number  #last lambda to try for given line, make sure to capture the whole line
        lambdas = absorption[np.logical_and(fill_in_here<absorption[:,0],absorption[:,0]<fill_in_here)] #array of lambda-values
        relevant_noise = noise[np.logical_and(fill_in_here<absorption[:,0],absorption[:,0]<fill_in_here)]  #noise over the given lambda-range
        min_value =  #smallest flux value to test
        flux_min = np.linspace(min_value,1,n_flux) #flux values to try
        
        #***********************************************************************
        #**********choose *one* of the two solutions:***************************
        #***********************************************************************

        #solution no.1, the elegant one:
        #***********************************************************************
        comb=np.array(np.meshgrid(flux_min,sigma,lambdas[:,0],)).T.reshape(-1,3) #meshgrid of all possible combinations of values
        #check out examples and explanations for use of meshgrid in the Numpy section of the Numerical Compendium.
        sums=np.zeros(len(comb))  #chi2-values to be stored here
        
        for i in range(len(comb)): #loop over all combinations
            sums[i]= #chi^2 value for a given combination

        j=(abs(sums)).argmin() #minimize chi2
        bestmodel = F_model(fill_in_here)
        
        print('velocity:',)
        print('temp:',)
        print('flux_min:',)
        
        # solution no.2, the easy-to-understand solution:
        #**********************************************************************
        smsmin = #initial value
        sgm_best = 0.
        flx_best = 0.
        lmb_best = 0.
        for flx in flux_min:
            for sgm in sigma:
                for lmb in lambdas[:,0]:
                    sms = #chi^2 value
                    if (find_a_logical_test_here):
                        #fill in

        bestmodel = F_model(fill_in_here)
        
        print('velocity:',)
        print('temp:',)
        print('flux_min:',)
        
        plt.plot(lambdas[:,0],lambdas[:,1], label="Meassured flux"); plt.legend()
        plt.plot(lambdas[:,0], bestmodel, label="Gaussian fit"); plt.legend()
        plt.legend(loc="best")
        plt.xlabel("nm")
        plt.ylabel("$W/m^2$")
        plt.title("%s with line %.i" %(Mn(lines),lines))
        plt.savefig("%s%.i" %(Mn(lines),lines))
        plt.show()

        
   
absorption_line()