#!/usr/bin/env python


# This contains routines
# that are usefull
# when running
# the Fortran version of ARTDECO

from __future__ import print_function, division, absolute_import
from builtins import range

import numpy as np
import pylab as pl
import string
import os

##############################################

wvl_rtol =  0.01 # %

def approx_equal(a, b, tol):
     return abs(a - b) < tol

##############################################


def get_betal_all(ptcle_type, filename):

    f   = open(filename, 'r')
    tmp = ' '

    while tmp.find(' number of wavelengths') == -1 :
        tmp = f.readline()
    tmp   = f.readline()
    nlamb = int(tmp)

    tmp   = f.readline() # skip a line
    lamb  = np.zeros(nlamb)
    cext  = np.zeros(nlamb)
    ssa   = np.zeros(nlamb)
    g     = np.zeros(nlamb)
    trunc_coef = np.zeros(nlamb)
    nbetal = np.zeros(nlamb, dtype=int)

    for i in range(nlamb):
        tmp = f.readline()
        lamb[i]  = float(tmp.split()[0])
        nbetal[i] = int(tmp.split()[1])
        trunc_coef[i] = float(tmp.split()[2])
        cext[i]  = float(tmp.split()[3])
        ssa[i]   = float(tmp.split()[4])
        g[i]     = float(tmp.split()[5])            

    nbetal_tmp = np.max(nbetal)
    alpha1 = np.zeros((nlamb,nbetal_tmp+1))
    alpha2 = np.zeros((nlamb,nbetal_tmp+1))
    alpha3 = np.zeros((nlamb,nbetal_tmp+1))
    alpha4 = np.zeros((nlamb,nbetal_tmp+1))
    beta1 = np.zeros((nlamb,nbetal_tmp+1))
    beta2 = np.zeros((nlamb,nbetal_tmp+1))
    
    tmp = f.readline()
    for i in range(nlamb): 
         tmp = f.readline()
         tmp = f.readline()
         for j in range(nbetal[i]+1):
             tmp = f.readline()
             alpha1[i,j]= float(tmp.split()[1])
             alpha2[i,j]= float(tmp.split()[2])
             alpha3[i,j]= float(tmp.split()[3])
             alpha4[i,j]= float(tmp.split()[4])
             beta1[i,j]= float(tmp.split()[5])
             beta2[i,j]= float(tmp.split()[6])
             
    f.close()

    return lamb, alpha1, alpha2, alpha3, alpha4, beta1, beta2, nbetal, trunc_coef   


################################################################################"

def get_opt_all(dir_opt, ptcle_type, recomp='null'):

    if recomp != 'null':
        opt_file    =  recomp
 
    else:
        opt_file    =  dir_opt+'/opt_' + ptcle_type +'.dat'

    f   = open(opt_file, 'r')
    tmp = ' '

    if recomp == 'null':
        while tmp.find(' phase matrix elements') == -1 :
            tmp = f.readline()
        tmp   = f.readline()
        nelem = int(tmp)
    else:
        nelem=6
 
    while tmp.find(' number of wavelengths') == -1 :
        tmp = f.readline()
    tmp   = f.readline()
    nlamb = int(tmp)
    tmp   = f.readline() # skip a line
    lamb  = np.zeros(nlamb)
    nteta = np.zeros(nlamb, dtype=int)
    cext  = np.zeros(nlamb)
    ssa   = np.zeros(nlamb)
    g     = np.zeros(nlamb)

    if recomp != 'null':
        sumf11 = np.zeros(nlamb)
        trunc_coef = np.zeros(nlamb)
        nbetal = np.zeros(nlamb, dtype=int)

    for i in range(nlamb):
        tmp = f.readline()
        lamb[i]  = float(tmp.split()[0])
        if recomp != 'null':
            nbetal[i] = int(tmp.split()[1])
            nteta[i]  = int(tmp.split()[2])
            trunc_coef[i] = float(tmp.split()[3])
            sumf11[i]= float(tmp.split()[4])
            cext[i]  = float(tmp.split()[5])
            ssa[i]   = float(tmp.split()[6])
            g[i]     = float(tmp.split()[7])            
        else: 
            nteta[i] = int(tmp.split()[1])
            cext[i]  = float(tmp.split()[2])
            ssa[i]   = float(tmp.split()[3])
            g[i]     = float(tmp.split()[4])


    if len(np.unique(nteta)) !=1 :
         print(" ( get_opt_all) error nteta must be all the same")
         exit()

    nang = nteta[0]
    u   = np.zeros((nlamb, nang))
    f11 = np.zeros((nlamb, nang))
    f44 = np.zeros((nlamb, nang))
    f21 = np.zeros((nlamb, nang))
    f34 = np.zeros((nlamb, nang))
    if(nelem==6):
        f22 = np.zeros((nlamb, nang))
        f33 = np.zeros((nlamb, nang))
    for iwvl in range(nlamb):
         skipcomment(f)
         for j in range(nteta[iwvl]):
             tmp = f.readline()
             u[iwvl,j]   = float(tmp.split()[0])
             if nelem==4:
                 f11[iwvl,j]= float(tmp.split()[1])
                 f44[iwvl,j]= float(tmp.split()[2])
                 f21[iwvl,j]= float(tmp.split()[3])
                 f34[iwvl,j]= float(tmp.split()[4])
             else:    
                 f11[iwvl,j]= float(tmp.split()[1])
                 f22[iwvl,j]= float(tmp.split()[2])
                 f33[iwvl,j]= float(tmp.split()[3])
                 f44[iwvl,j]= float(tmp.split()[4])
                 f21[iwvl,j]= float(tmp.split()[5])
                 f34[iwvl,j]= float(tmp.split()[6])
    f.close()

    for j in range(nang):
         if len(np.unique(u[:,j])) !=1 :
              print(" ( get_opt_all) error teta must be all the same")
              exit()

    u = u[0,:]     
    ang = np.arccos(u) / np.pi * 180.0

    if recomp != 'null': 
         return lamb, ang, f11, f22, f33, f44, f21, f34, cext, ssa, g, nbetal, nteta, trunc_coef, sumf11
    else:
        if nelem==4:
            return lamb, ang, f11, f44, f21, f34, cext, ssa, g
        else:
            return lamb, ang, f11, f22, f33, f44, f21, f34, cext, ssa, g

###############################################################################"

def get_opt(dir_opt, ptcle_type, wavel, recomp='null'):
 
    wvl_tol =  wavel *  ( wvl_rtol / 100.0)
    
    #print  wavel, wavel+wvl_tol
    #print "get_opt..."          

    flag_wl = False

    if recomp != 'null':
        opt_file    =  recomp
 
    else:
        opt_file    =  dir_opt+'/opt_' + ptcle_type +'.dat'

    f   = open(opt_file, 'r')
    tmp = ' '

    if recomp == 'null':
        while tmp.find(' phase matrix elements') == -1 :
            tmp = f.readline()
        tmp   = f.readline()
        nelem = int(tmp)
    else:
        nelem=6
 
    while tmp.find(' number of wavelengths') == -1 :
        tmp = f.readline()
    tmp   = f.readline()
    nlamb = int(tmp)
    tmp   = f.readline() # skip a line
    lamb  = np.zeros(nlamb)
    nteta = np.zeros(nlamb, dtype=int)
    cext  = np.zeros(nlamb)
    ssa   = np.zeros(nlamb)
    g     = np.zeros(nlamb)

    if recomp != 'null':
        sumf11 = np.zeros(nlamb)
        trunc_coef = np.zeros(nlamb)
        nbetal = np.zeros(nlamb, dtype=int)

    for i in range(nlamb):
        tmp = f.readline()
        lamb[i]  = float(tmp.split()[0])
        if approx_equal(lamb[i], wavel, wvl_tol):
            flag_wl = True
            wlindex = i
        if recomp != 'null':
            nbetal[i] = int(tmp.split()[1])
            nteta[i]  = int(tmp.split()[2])
            trunc_coef[i] = float(tmp.split()[3])
            sumf11[i]= float(tmp.split()[4])
            cext[i]  = float(tmp.split()[5])
            ssa[i]   = float(tmp.split()[6])
            g[i]     = float(tmp.split()[7])            
        else: 
            nteta[i] = int(tmp.split()[1])
            cext[i]  = float(tmp.split()[2])
            ssa[i]   = float(tmp.split()[3])
            g[i]     = float(tmp.split()[4])
    if (flag_wl==False):
        print(' Error :')
        print('    The wavelength '+'%.5f'%wavel+' microns is not stored in opt_file')
        print('')
        exit()
    nang = nteta[wlindex]
    u   = np.zeros(nang)
    f11 = np.zeros(nang)
    f44 = np.zeros(nang)
    f21 = np.zeros(nang)
    f34 = np.zeros(nang)
    if(nelem==6):
        f22 = np.zeros(nang)
        f33 = np.zeros(nang)
    for i in range(wlindex):
        skipcomment(f)
        for j in range(nteta[i]):
            tmp = f.readline()
    skipcomment(f)
    for j in range(nteta[wlindex]):
        tmp = f.readline()
        u[j]   = float(tmp.split()[0])
        if nelem==4:
            f11[j]= float(tmp.split()[1])
            f44[j]= float(tmp.split()[2])
            f21[j]= float(tmp.split()[3])
            f34[j]= float(tmp.split()[4])
        else:    
            f11[j]= float(tmp.split()[1])
            f22[j]= float(tmp.split()[2])
            f33[j]= float(tmp.split()[3])
            f44[j]= float(tmp.split()[4])
            f21[j]= float(tmp.split()[5])
            f34[j]= float(tmp.split()[6])
    f.close()
    ang = np.arccos(u) / np.pi * 180.0

    if recomp != 'null': 
        return ang, f11, f22, f33, f44, f21, f34, cext[wlindex], ssa[wlindex], g[wlindex], nbetal[wlindex], nteta[wlindex], trunc_coef[wlindex], sumf11[wlindex]                      
    else:
        if nelem==4:
            return ang, f11, f44, f21, f34, cext[wlindex], ssa[wlindex], g[wlindex]
        else:
            return ang, f11, f22, f33, f44, f21, f34, cext[wlindex], ssa[wlindex], g[wlindex]

##############################################

def get_betal(ptcle_type, wavel, filename):
    
    wvl_tol =  wavel *  ( wvl_rtol / 100.0)
    
    flag_wl = False

    f   = open(filename, 'r')
    tmp = ' '

    while tmp.find(' number of wavelengths') == -1 :
        tmp = f.readline()
    tmp   = f.readline()
    nlamb = int(tmp)

    tmp   = f.readline() # skip a line
    lamb  = np.zeros(nlamb)
    cext  = np.zeros(nlamb)
    ssa   = np.zeros(nlamb)
    g     = np.zeros(nlamb)
    trunc_coef = np.zeros(nlamb)
    nbetal = np.zeros(nlamb, dtype=int)

    for i in range(nlamb):
        tmp = f.readline()
        lamb[i]  = float(tmp.split()[0])
        if approx_equal(lamb[i], wavel, wvl_tol):
            flag_wl = True
            wlindex = i
        nbetal[i] = int(tmp.split()[1])
        trunc_coef[i] = float(tmp.split()[2])
        cext[i]  = float(tmp.split()[3])
        ssa[i]   = float(tmp.split()[4])
        g[i]     = float(tmp.split()[5])            

        #print lamb[i], nbetal[i] 

    if (flag_wl==False):
        print(' Error :')
        print('    The wavelength '+'%.5f'%wavel+' microns is not stored in opt_file')
        print('')
        exit()
    nbetal_tmp = nbetal[wlindex]
    alpha1 = np.zeros(nbetal_tmp+1)
    alpha2 = np.zeros(nbetal_tmp+1)
    alpha3 = np.zeros(nbetal_tmp+1)
    alpha4 = np.zeros(nbetal_tmp+1)
    beta1 = np.zeros(nbetal_tmp+1)
    beta2 = np.zeros(nbetal_tmp+1)
    
    # skip lines
    tmp = f.readline()
    #print tmp
    for i in range(wlindex):
        tmp = f.readline()
        tmp = f.readline()
        for j in range(nbetal[i]+1):
            tmp = f.readline()
            #print j, tmp

    #print ''        
    #print wavel    
    tmp = f.readline()
    #print tmp
    #print ''

    tmp = f.readline()
    for j in range(nbetal[wlindex]+1):
        tmp = f.readline()
        alpha1[j]= float(tmp.split()[1])
        alpha2[j]= float(tmp.split()[2])
        alpha3[j]= float(tmp.split()[3])
        alpha4[j]= float(tmp.split()[4])
        beta1[j]= float(tmp.split()[5])
        beta2[j]= float(tmp.split()[6])
    f.close()

    return alpha1, alpha2, alpha3, alpha4, beta1, beta2, nbetal[wlindex], trunc_coef[wlindex]    

##############################################

def set_nstr(rt_model, nstr, dir_artdeco):
    dir_artdeco = dir_artdeco+'input/' 
    if rt_model == 'doad':
        file = dir_artdeco+'doad_spec.dat'
    elif rt_model == 'addoub':
        file = dir_artdeco+'ad_spec.dat'
    elif rt_model == 'disort':
        file = dir_artdeco+'od_spec.dat'
    else:
        print('(set_nstr) Error :')
        print('(set_nstr) Model name not recognized')
        return
    # read file
    f = open(file, 'r')
    file_content = f.readlines()
    f.close
    # change nstr
    i=0
    while file_content[i].find('nstr') == -1 :
        i=i+1
    file_content[i+1] = str(int(nstr))+'  \n'
    # write result
    f = open(file, 'w')
    for i in range(len(file_content)):
        f.write(file_content[i])
    f.close
    

####################################################

def read_solrad_out(res_file):

    tmp = ''
    f = open(res_file, 'r')
    while tmp.find('Number of wavelengths') == -1 :
        tmp = f.readline()
    tmp  = f.readline()
    nwvl =  int(tmp)
    wvl  = np.zeros(nwvl) 
    F0   = np.zeros(nwvl)
    tmp  = f.readline() # comment
    tmp  = f.readline()
    nsza =  int(tmp)
    sza    = np.zeros(nsza) 
    varsol = np.zeros(nsza)
    lat    = np.zeros(nsza)
    lon    = np.zeros(nsza)
    doy    = np.zeros(nsza, dtype=int)
    h_tu   = np.zeros(nsza)
    tmp  = f.readline() # comment
    tmp  = f.readline() # comment
    for i in range(nwvl):
        tmp  = f.readline()
        #print tmp
        wvl[i] = float(tmp.split()[0])
        F0[i] = float(tmp.split()[1])
    tmp  = f.readline() # comment
    tmp  = f.readline() # comment
    for i in range(nsza):
        tmp  = f.readline()
        #print tmp
        sza[i]    = float(tmp.split()[0])
        varsol[i] = float(tmp.split()[1])
        lon[i]    = float(tmp.split()[2])
        lat[i]    = float(tmp.split()[3])
        doy[i]    = int(tmp.split()[4])
        h_tu[i]   = float(tmp.split()[5])
    f.close()

    return [wvl,F0,sza,varsol,lon,lat,doy,h_tu]
    
    
####################################################

def read_flux_kdis_band(res_file):

    tmp =  ' '
    
    f = open(res_file, 'r')

    while tmp.find('wavelengths') == -1 :
        tmp = f.readline()
    tmp   = f.readline()
    nwvl = int(tmp)
    #print nwvl

    while tmp.find('zenith') == -1 :
        tmp = f.readline()
    tmp   = f.readline()
    nsza = int(tmp)
    #print nsza

    while tmp.find('altitude') == -1 :
        tmp = f.readline()
    tmp   = f.readline()
    nalt = int(tmp)
    #print nalt
    tmp   = f.readline()

    wvl  = np.zeros(nwvl)
    sza  = np.zeros(nsza)
    alt  = np.zeros(nalt)
    flux = np.zeros((nsza, nwvl, nalt, 3))

    for iwvl in range(nwvl):
        tmp = f.readline()
        tmp = f.readline()
        wvl[iwvl] = float(tmp.split()[0])
        for isza in range(nsza):
            tmp = f.readline()
            tmp = f.readline()    
            if iwvl == 0 :
                sza[isza] =  float(tmp.split()[0])
            tmp = f.readline()    
            tmp = f.readline()    
            for ialt in range(nalt):
                tmp = f.readline()    
                if (ialt==0 and iwvl==0):
                    alt[ialt] = float(tmp.split()[0])
                flux[isza,iwvl,ialt,0] = float(tmp.split()[2])
                flux[isza,iwvl,ialt,1] = float(tmp.split()[3])
                flux[isza,iwvl,ialt,2] = float(tmp.split()[4])

    return[nwvl, nsza, nalt, wvl, sza, alt, flux]


####################################################

def read_artdeco_rad(res_file):

    tmp =  ' '
    
    f = open(res_file, 'r')
    while tmp.find('# total number of geometries') == -1 :
        tmp = f.readline()
    tmp   = f.readline()
    ngeom = int(tmp)
    while tmp.find('# nsza, nvza, nvaa') == -1 :
        tmp = f.readline()
    tmp  = f.readline()
    nsza = int(tmp.split()[0])
    nvza = int(tmp.split()[1])
    nvaa = int(tmp.split()[2])
    while tmp.find('# 4 = full polarization, 3 = 3x3 approximation, 1 = scalar') == -1 :
        tmp = f.readline()
    tmp  = f.readline()
    nmat = int(tmp)
    while tmp.find('# number of wavelengths') == -1 :
        tmp = f.readline()
    tmp  = f.readline()
    nwvl = int(tmp)
    
    if nmat == 1:
        SvR  = np.zeros((nsza, nvza, nvaa, nwvl))
        dSvR = np.zeros((nsza, nvza, nvaa, nwvl))
        polar = 0
    else:
        SvR  = np.zeros((nsza, nvza, nvaa, nmat, nwvl))
        dSvR = np.zeros((nsza, nvza, nvaa, nmat, nwvl))
        polar  = np.zeros((nsza, nvza, nvaa, nwvl))
    sza = np.zeros(nsza)
    vza = np.zeros(nvza)
    vaa = np.zeros(nvaa)
    wvl = np.zeros(nwvl)
    wvl     = np.zeros(nwvl)
    cputime = np.zeros(nwvl)
    theta = np.zeros((nsza, nvza, nvaa))

    tmp = f.readline() # skip comment line
    for iwvl in range(nwvl):
        tmp = f.readline() # skip comment line
        tmp = f.readline()
        wvl[iwvl] = float(tmp.split()[0])        
        tmp = f.readline() # skip comment line
        tmp = f.readline()
        cputime[iwvl] = float(tmp.split()[0])
        tmp = f.readline() # skip comment line
        for i in range(nsza):
            for j in range(nvza):
                for k in range(nvaa):
                    tmp = f.readline()
                    sza[i] = float(tmp.split()[0])
                    vza[j] = float(tmp.split()[1])
                    vaa[k] = float(tmp.split()[2])
                    theta[i,j,k] = float(tmp.split()[3])
                    if nmat == 1 :
                        SvR[i,j,k,iwvl]    = float(tmp.split()[4])
                        dSvR[i,j,k,iwvl]   = float(tmp.split()[5])
                    if nmat == 3 :
                        SvR[i,j,k,0,iwvl]  = float(tmp.split()[4])
                        SvR[i,j,k,1,iwvl]  = float(tmp.split()[5])
                        SvR[i,j,k,2,iwvl]  = float(tmp.split()[6])
                        dSvR[i,j,k,0,iwvl] = float(tmp.split()[7])
                        dSvR[i,j,k,1,iwvl] = float(tmp.split()[8])
                        dSvR[i,j,k,2,iwvl] = float(tmp.split()[9])
                        polar[i,j,k,iwvl]  = float(tmp.split()[10])
                    if nmat == 4:
                        SvR[i,j,k,0,iwvl]  = float(tmp.split()[4])
                        SvR[i,j,k,1,iwvl]  = float(tmp.split()[5])
                        SvR[i,j,k,2,iwvl]  = float(tmp.split()[6])
                        SvR[i,j,k,3,iwvl]  = float(tmp.split()[7])
                        dSvR[i,j,k,0,iwvl] = float(tmp.split()[8])
                        dSvR[i,j,k,1,iwvl] = float(tmp.split()[9])
                        dSvR[i,j,k,2,iwvl] = float(tmp.split()[10])
                        dSvR[i,j,k,3,iwvl] = float(tmp.split()[11])
                        polar[i,j,k,iwvl]  = float(tmp.split()[12])
    f.close()
   
    return [nmat, sza, vza, vaa, wvl, theta, SvR, dSvR, polar, cputime]

####################################################

def read_artdeco_rad_filter(res_files):

    tmp =  ' '

    nfilter = len(res_files)

    for ifilt in range(nfilter) :
        
        f = open(res_files[ifilt], 'r')
        while tmp.find('# total number of geometries') == -1 :
            tmp = f.readline()
        tmp   = f.readline()
        ngeom = int(tmp)
        while tmp.find('# nsza, nvza, nvaa') == -1 :
            tmp = f.readline()
        tmp  = f.readline()
        nsza = int(tmp.split()[0])
        nvza = int(tmp.split()[1])
        nvaa = int(tmp.split()[2])
        while tmp.find('# 4 = full polarization, 3 = 3x3 approximation, 1 = scalar') == -1 :
            tmp = f.readline()
        tmp  = f.readline()
        nmat = int(tmp)

        if (ifilt == 0) :
            if nmat == 1:
                SvR  = np.zeros((nsza, nvza, nvaa, nfilter))
                dSvR = np.zeros((nsza, nvza, nvaa, nfilter))
                polar = 0
            else:
                SvR  = np.zeros((nsza, nvza, nvaa, nmat, nfilter))
                dSvR = np.zeros((nsza, nvza, nvaa, nmat, nfilter))
                polar  = np.zeros((nsza, nvza, nvaa, nfilter))
            sza = np.zeros(nsza)
            vza = np.zeros(nvza)
            vaa = np.zeros(nvaa)
            cputime = np.zeros(nfilter)
            theta = np.zeros((nsza, nvza, nvaa))

        tmp = f.readline() # skip comment line
        tmp = f.readline()
        cputime[ifilt] = float(tmp.split()[0])
        tmp = f.readline()
        for i in range(nsza):
            for j in range(nvza):
                for k in range(nvaa):
                    tmp = f.readline()
                    sza[i] = float(tmp.split()[0])
                    vza[j] = float(tmp.split()[1])
                    vaa[k] = float(tmp.split()[2])
                    theta[i,j,k] = float(tmp.split()[3])
                    if nmat == 1 :
                        SvR[i,j,k,ifilt]    = float(tmp.split()[4])
                        dSvR[i,j,k,ifilt]   = float(tmp.split()[5])
                    if nmat == 3 :
                        SvR[i,j,k,0,ifilt]  = float(tmp.split()[4])
                        SvR[i,j,k,1,ifilt]  = float(tmp.split()[5])
                        SvR[i,j,k,2,ifilt]  = float(tmp.split()[6])
                        dSvR[i,j,k,0,ifilt] = float(tmp.split()[7])
                        dSvR[i,j,k,1,ifilt] = float(tmp.split()[8])
                        dSvR[i,j,k,2,ifilt] = float(tmp.split()[9])
                        polar[i,j,k,ifilt]  = float(tmp.split()[10])
                    if nmat == 4:
                        SvR[i,j,k,0,ifilt]  = float(tmp.split()[4])
                        SvR[i,j,k,1,ifilt]  = float(tmp.split()[5])
                        SvR[i,j,k,2,ifilt]  = float(tmp.split()[6])
                        SvR[i,j,k,3,ifilt]  = float(tmp.split()[7])
                        dSvR[i,j,k,0,ifilt] = float(tmp.split()[8])
                        dSvR[i,j,k,1,ifilt] = float(tmp.split()[9])
                        dSvR[i,j,k,2,ifilt] = float(tmp.split()[10])
                        dSvR[i,j,k,3,ifilt] = float(tmp.split()[11])
                        polar[i,j,k,ifilt]  = float(tmp.split()[12])
        f.close()
    
    return [nmat, sza, vza, vaa, theta, SvR, dSvR, polar, cputime]

##############################################

def set_artdeco_in(keywords, saveroot, mode, filters, wavel,           \
                       ptcle, trunc_model, nbetal, rt_model,           \
                       corint, thermal, nmat, atm, surface, surf_temp, \
                       solar_spec, sza, vza, vaa):

    # This routine returns the artdeco_in.dat file 
    # content

    nfilters = len(filters)
    nptcle   = len(ptcle)
    nsza     = len(sza)
    nvza     = len(vza)
    nvaa     = len(vaa)

    artdeco_in=['# \n']
    artdeco_in.append('# MAIN INPUT FILE FOR ARTDECO PROGRAM  \n')
    artdeco_in.append('#   \n')
    artdeco_in.append('######################## \n')
    artdeco_in.append('# keywords (Ex: none, verbose,...) \n')
    artdeco_in.append(keywords+' \n')
    artdeco_in.append('######################## \n')
    artdeco_in.append('# outfiles root name \n')
    artdeco_in.append(saveroot+' \n')
    artdeco_in.append('######################## \n')
    artdeco_in.append('# mode \n')
    artdeco_in.append(mode+' \n')
    artdeco_in.append('#######################\n')
    artdeco_in.append('#     filters\n')
    if nfilters==0:
         artdeco_in.append('  none  \n')
    else :
         for ifilt in range(nfilters):
              artdeco_in.append(filters[ifilt]+' \n')
    artdeco_in.append('######################## \n')
    artdeco_in.append('# Wavelengths (microns) \n')
    s= ''
    for iwl in range(len(wavel)):
        s = s+'  %.5f'%wavel[iwl]
    artdeco_in.append(s+'   \n') 
    artdeco_in.append('###################### \n')
    artdeco_in.append('# Particles \n')
    artdeco_in.append('# type           interp.       H-G         Tau_550          vertical distribution type     vdist parameters (km)  \n')
    if nptcle == 0 :
         artdeco_in.append(' none \n')
    else :
         for iptcle in range(nptcle): 
              artdeco_in.append(ptcle[iptcle]+'   \n')
    artdeco_in.append('# \n')
    artdeco_in.append('########################## \n')
    artdeco_in.append('# Truncation method (none, dfit, DM, Potter) \n')
    artdeco_in.append(trunc_model+'    \n')
    artdeco_in.append('########################## \n')
    artdeco_in.append('# Number of Betal \n')
    artdeco_in.append('%4i'%nbetal+'   \n')
    artdeco_in.append('########################## \n')
    artdeco_in.append('# Radiative transfer model \n')
    artdeco_in.append(rt_model+'   \n')
    artdeco_in.append('#########################  \n')
    artdeco_in.append('# TMS correction (yes/no)   \n')
    artdeco_in.append(corint+'    \n')
    artdeco_in.append('######################### \n')
    artdeco_in.append('# Thermal (yes/no) \n')
    artdeco_in.append(thermal+' \n')
    artdeco_in.append('######################### \n')
    artdeco_in.append('# nmat: polarization included (4) or scalar (1) or (3) for I, Q, U \n')
    artdeco_in.append('%4i'%nmat+' \n')
    artdeco_in.append('######################### \n')
    artdeco_in.append('# Atmosphere model  \n')
    artdeco_in.append(atm+' \n')
    artdeco_in.append('######################## \n')
    artdeco_in.append('# Surface \n')
    artdeco_in.append('# family     type     albedo (if family=lambert & type=cste) \n')
    artdeco_in.append(surface+' \n')
    artdeco_in.append('# Surface temperature (in K) (-1: to use the same temperature as the atmosphere bottom temperature)\n')
    artdeco_in.append('%.4f'%surf_temp+' \n')
    artdeco_in.append('######################## \n')
    artdeco_in.append('#  Incident radiation characterisrtics \n')
    artdeco_in.append('# spectrum :\n')
    artdeco_in.append(solar_spec +' \n')
    artdeco_in.append('# solar zenith angle(s) : between 0.0 and 90.0 degrees \n')
    s= ''
    for isza in range(nsza):
        s = s+'  %.4f'%sza[isza]
    artdeco_in.append(s+'   \n')     
    artdeco_in.append('# Longitude, Latitude, day of the year,  time (decimal, TU) \n')
    artdeco_in.append(' none \n')
    artdeco_in.append('########################### \n')
    artdeco_in.append('# Geometrical conditions : \n')
    artdeco_in.append('# view zenith angle(s) : between 90.0 (exluded) and 0.0 degrees \n')
    s= ''
    for ivza in range(nvza):
        s = s+'  %.4f'%vza[ivza]
    artdeco_in.append(s+'   \n') 
    artdeco_in.append('# view azimuthal angle(s) : between 0. and 180. degrees (0 : specular direction, 180 : backscattering direction) \n')
    s= ''
    for ivaa in range(nvaa):
        s = s+'  %.4f'%vaa[ivaa]
    artdeco_in.append(s+'   \n')

    return artdeco_in


########################################################################################

def skipcomment(f):
    pos = f.tell()
    tmp = f.readline()
    while tmp.startswith('#'):
        pos = f.tell()
        tmp = f.readline()
    f.seek(pos)


if __name__=='__main__':
    print("pyartdeco_utils is a library")