AiresInfoFunctionsGRANDROOT.py

##
## WARNING WARNING WARNING
##
#this file is a local copy of AiresInfoFunctions from https://github.com/mjtueros/ZHAireS-Python  (on Jun 25 2021)so that 
#you dont have to setup that repository, that requires ZHAireS to be installed. (since we dont use it for now). Will update it manually if needed
#table exporting will not be available!!
##
## WARNING WARNING WARNING
##


#if you are having encoding errors, specially in lyon, be sure to force the environment to utf8
#setenv  LANGUAGE en_US.UTF-8 ;setenv  LC_ALL en_US.UTF-8 ;setenv  LANG en_US.UTF-8setenv LC_TYPE en_US.UTF-8
#

#this functions will accept GRAND and AIRES outmode, to give the results in each convention
#it will output the primary zen,azim,energy,primarytype, taken from the .inp file present at input_file_path) (assumed only one .inp file per dir)


#TO DO: treat correctly the different possible primary types, including RASPASS Multi primary
#TO DO: treat correctly the case where a distribution of primary types or energies or angles is put in the input
#TO DO: Check consistency of the output (energy within a range, angles within a range, etc)
#TO DO: have the .sry reader regenerate the summary using AiresSry, if the file is not foun

#6/2019 Matias Tueros, first attempt at python based on original script by Anne Zilles.
#10/2019 Migrated them to make it the official library
#12/2019 Set up on git
#03/2020 Corona Virus

import sys
from sys import argv
import os

import glob
import logging

#this function reads he Aires .sry file and tries to extract information from the simulation parameters
#using the .sry file is preferred to using the .inp files because:
#Output is standarized: you dont know what you can find in an .inp file (leading spaces, commented lines, repeated keywords, etc)
#Output is what really happened, not what the user whished it would happen when he did his crappy .inp file.

#lets break the ReadAiresSry into smaller modules. It has the disadventage of opening, scanning and closing the file each time
#but it adds modularity, and closes the files when it does not use it any more.
#If at some point we need speed, then we could input datafile, and make a wraper for opening the file

AiresPath="/home/mjtueros/aires/bin"
#AiresPath=os.environ["AIRESBINDIR"]
#AiresPath="/home/mjtueros/AiresRepository/Dropbox/AiresBzr/Aires.19-04-04-ZHAireS-development/zhadl/bin"

def GetZenithAngleFromSry(sry_file,outmode="GRAND"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Primary zenith angle:' in line:
          line = line.lstrip()
          stripedline=line.split(' ',-1)
          zen=float(stripedline[len(stripedline)-2])
          if outmode == 'GRAND':
            zen = 180-zen  #conversion to GRAND convention i.e. pointing towards antenna/propagtion direction
          #logging.debug('Found Zenith ' + str(zen))
          return zen
      try:
        zen
      except NameError:
        zen = 0 #If no zenith angle was included in the input file, AIRES defaults to 0
        if outmode == 'GRAND':
          zen = 180-0 #that translates to 180 in GRAND
        logging.info("Zenith Angle not found in sry file, defaulting to:" + str(zen))
        return zen
  except:
    logging.error("GetZenithAngleFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetAzimuthAngleFromSry(sry_file,outmode="GRAND"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Primary azimuth angle:' in line:
          line = line.lstrip()
          stripedline=line.split(' ',-1)
          azim = float(stripedline[len(stripedline)-2])
          if outmode == 'GRAND':
            azim=azim+180 #conversion to GRAND convention i.e. pointing towards antenna/propagtion direction
            if azim>=360:
              azim= azim-360
          #logging.debug('Found Azimuth ' + str(azim))
          return azim
      try:
        azim
      except NameError:
        azim = 0 #If no azimuth angle was included in the input file, AIRES defaults to 0
        if outmode == 'GRAND':
          azim = 0+180 # that translates to 18
        logging.info("Azimuth Angle not found in sry file, defaulting to:" + str(azim))
        return azim

  except:
    logging.error("GetAzimuthAngleFromSry:file not found or invalid:"+sry_file)
    raise
    return -1


def GetEnergyFromSry(sry_file,outmode="GRAND"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Primary energy:' in line:
          line = line.lstrip()
          stripedline=line.split(' ',-1)
          try: #this is to detect if it is the first time we find the string, so that we ignore the following
            energy
          except NameError:
            energy = float(stripedline[len(stripedline)-2])
            unit= str(stripedline[len(stripedline)-1])
            if unit == "eV\n":
             energy = energy *1e-9
            if unit == "KeV\n":
             energy = energy *1e-6
            if unit == "MeV\n":
              energy = energy *1e-3
            if unit == "GeV\n":
              energy = energy
            if unit == "TeV\n":
              energy = energy *1e3
            if unit == "PeV\n":
              energy = energy *1e6
            if unit == "EeV\n":
              energy = energy *1e9

            if outmode == 'GRAND': #AIRES mode outputs in GeV, GRAND in EeV
              energy = energy * 1e-9
            #logging.debug('Found Energy ' + str(energy)) #debug level 1
            return energy
      try:
        energy
      except NameError:
        logging.error('warning energy not found, Aires has no default value,  cannot continue')
        exit()
  except:
    logging.error("GetEnergyFromSry:file not found or invalid:"+sry_file)
    raise
    return -1


#output is in meters
def GetCorePositionFromInp(inp_file,outmode="N/A"):
  try:
    datafile=open(inp_file,'r')
    with open(inp_file, "r") as datafile:
      for line in datafile:
        if '#Core Position:' in line:
          line = line.lstrip()
          stripedline=line.split(':',-1)
          stripedline=stripedline[1]
          stripedline=stripedline.split(' ',-1)
          x=float(stripedline[1])
          y=float(stripedline[2])
          z=float(stripedline[3])
          coreposition=(x,y,z)
          return coreposition
      try:
        coreposition
      except NameError:
        logging.error('warning core position not found, defaulting to (0,0,0)')
        return (0.0,0.0,0.0)
  except:
    logging.error("GetCorePositionFromInp:file not found or invalid:"+inp_file)
    raise
    return -1



def GetThinningRelativeEnergyFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Thinning energy:' in line:
          line = line.lstrip()
          stripedline=line.split('Thinning energy:',-1)
          stripedline=stripedline[1]
          stripedline=stripedline.lstrip()
          stripedline=stripedline.split(' ',-1)
          Thinning=stripedline[0]
          Relative=stripedline[1]
          if(Relative!="Relative\n"):
             logging.error('warning , we only support Relative thinning for now, sorry!. This is easy to implement if you need it!')
             return -1
          return Thinning
      try:
        Thinning
      except NameError:
        logging.error('warning Thinning energy not found, Aires has no default value, cannot continue')
        exit()
  except:
    logging.error("GetThinningRelativeEnergyFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetGammaEnergyCutFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Cut energy for gammas:' in line:
          line = line.lstrip()
          stripedline=line.split('Cut energy for gammas:',-1)
          stripedline=stripedline[1]
          stripedline=stripedline.lstrip()
          stripedline=stripedline.split(' ',-1)
          energy=stripedline[0]
          unit=stripedline[1]
          if unit == "eV\n":
            energy = energy *1e-6
          if unit == "KeV\n":
            energy = energy *1e-3
          if unit == "MeV\n":
            energy = energy
          if unit == "GeV\n":
            energy = energy *1e3
          if unit == "TeV\n":
            energy = energy *1e6
          if unit == "PeV\n":
            energy = energy *1e9
          if unit == "EeV\n":
            energy = energy *1e12
          return energy
      try:
        energy
      except NameError:
        logging.error('warning GammaCutEnergy not found, defaulting to 80keV')
        return 80+1e-3
  except:
    logging.error("GetGammaEnergyCutFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetElectronEnergyCutFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Cut energy for e+ e-:' in line:
          line = line.lstrip()
          stripedline=line.split('Cut energy for e+ e-:',-1)
          stripedline=stripedline[1]
          stripedline=stripedline.lstrip()
          stripedline=stripedline.split(' ',-1)
          energy=stripedline[0]
          unit=stripedline[1]
          if unit == "eV\n":
            energy = energy *1e-6
          if unit == "KeV\n":
            energy = energy *1e-3
          if unit == "MeV\n":
            energy = energy
          if unit == "GeV\n":
            energy = energy *1e3
          if unit == "TeV\n":
            energy = energy *1e6
          if unit == "PeV\n":
            energy = energy *1e9
          if unit == "EeV\n":
            energy = energy *1e12
          return energy
      try:
        energy
      except NameError:
        logging.error('warning ElectronCutEnergy not found, defaulting to 80keV')
        return 80+1e-3
  except:
    logging.error("GetElectronEnergyCutFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetMuonEnergyCutFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Cut energy for mu+ mu-:' in line:
          line = line.lstrip()
          stripedline=line.split('Cut energy for mu+ mu-:',-1)
          stripedline=stripedline[1]
          stripedline=stripedline.lstrip()
          stripedline=stripedline.split(' ',-1)
          energy=stripedline[0]
          unit=stripedline[1]
          if unit == "eV\n":
            energy = energy *1e-6
          if unit == "KeV\n":
            energy = energy *1e-3
          if unit == "MeV\n":
            energy = energy
          if unit == "GeV\n":
            energy = energy *1e3
          if unit == "TeV\n":
            energy = energy *1e6
          if unit == "PeV\n":
            energy = energy *1e9
          if unit == "EeV\n":
            energy = energy *1e12
          return energy
      try:
        energy
      except NameError:
        logging.error('warning MuonCutEnergy not found, defaulting to 10MeV')
        return 10
  except:
    logging.error("GetMuonEnergyCutFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetMesonEnergyCutFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Cut energy for mesons:' in line:
          line = line.lstrip()
          stripedline=line.split('Cut energy for mesons:',-1)
          stripedline=stripedline[1]
          stripedline=stripedline.lstrip()
          stripedline=stripedline.split(' ',-1)
          energy=stripedline[0]
          unit=stripedline[1]
          if unit == "eV\n":
            energy = energy *1e-6
          if unit == "KeV\n":
            energy = energy *1e-3
          if unit == "MeV\n":
            energy = energy
          if unit == "GeV\n":
            energy = energy *1e3
          if unit == "TeV\n":
            energy = energy *1e6
          if unit == "PeV\n":
            energy = energy *1e9
          if unit == "EeV\n":
            energy = energy *1e12
          return energy
      try:
        energy
      except NameError:
        logging.error('warning MesonCutEnergy not found, defaulting to 60MeV')
        return 60
  except:
    logging.error("GetMesonEnergyCutFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetNucleonEnergyCutFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Cut energy for nucleons:' in line:
          line = line.lstrip()
          stripedline=line.split('Cut energy for nucleons:',-1)
          stripedline=stripedline[1]
          stripedline=stripedline.lstrip()
          stripedline=stripedline.split(' ',-1)
          energy=stripedline[0]
          unit=stripedline[1]
          if unit == "eV\n":
            energy = energy *1e-6
          if unit == "KeV\n":
            energy = energy *1e-3
          if unit == "MeV\n":
            energy = energy
          if unit == "GeV\n":
            energy = energy *1e3
          if unit == "TeV\n":
            energy = energy *1e6
          if unit == "PeV\n":
            energy = energy *1e9
          if unit == "EeV\n":
            energy = energy *1e12
          return energy
      try:
        energy
      except NameError:
        logging.error('warning NucleonCutEnergy not found, defaulting to 60MeV')
        return 60
  except:
    logging.error("GetNucleonEnergyCutFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetPrimaryFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Primary particle:' in line:
          line = line.lstrip()
          stripedline=line.split(' ',-1)
          if(len(stripedline)==3):
            primarytype = str(stripedline[len(stripedline)-1])
          elif(len(stripedline)==5):
            primarytype = str(stripedline[len(stripedline)-3])
          elif(len(stripedline)==6):
            primarytype = str(stripedline[len(stripedline)-4])
          elif(len(stripedline)==7):
            primarytype = str(stripedline[len(stripedline)-5])
          else:
            primarytype = "unknown"
          primarytype=primarytype.replace('\n','')
          #logging.debug('Found Primary ' + primarytype) #debug level 1

          if(outmode=="GRAND"):
           #convert to GRAND coding. PDG?
            if primarytype=="Proton":
              primarytype=2212
            if primarytype=="Neutron":
              primarytype=2112
            if primarytype=="Pi0":
              primarytype=111
            if primarytype=="Pi+":
              primarytype=211
            if primarytype=="Pi-":
              primarytype=-211

          return primarytype
      try:
        primarytype
      except NameError:
        logging.error('warning primary not found, Aires has no default value, cannot continue')
        exit()
  except:
    logging.error("GetPrimaryFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetSlantXmaxFromSry(sry_file,outmode="N/A"): #To do. Handle when Xmax is not found, becouse the fit didnt converge
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Sl. depth of max. (g/cm2)' in line:
          line = line.lstrip()
          stripedline=line.split(' ',-1)
          xmax = float(stripedline[len(stripedline)-1])
          #logging.debug('Found Xmax ' + str(xmax)) #debug level 1
          return xmax
      try:
        xmax
      except NameError:
        logging.info('warning xmax not found')
        xmax=-1
        return xmax
  except:
    logging.error("GetSlantXmaxFromSry:file not found or invalid:"+sry_file)
    raise
    return -1
#
#Charged pcles. at maximum
def GetNmaxFromSry(sry_file,outmode="N/A"): #To do. Handle when Xmax is not found, becouse the fit didnt converge
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Charged pcles. at maximum' in line:
          line = line.lstrip()
          stripedline=line.split(' ',-1)
          nmax = float(stripedline[len(stripedline)-1])
          #logging.debug('Found Nmax ' + str(nmax)) #debug level 1
          return nmax
      try:
        xmax
      except NameError:
        logging.info('warning nmax not found')
        xmax=-1
        return xmax
  except:
    logging.error("GetNmaxFromSry:file not found or invalid:"+sry_file)
    raise
    return -1
#                              Altitude  Distance     x        y        z
#      Location of max.(Km):     3.612     3.66     0.00     0.57     3.61
def GetKmXmaxFromSry(sry_file,outmode="N/A"): #To do. Handle when Xmax is not found, becouse the fit didnt converge, or becouse this is not ZHAireS, or its latest version
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if ('Location of max.(Km)' in line) or ('Pos. Max.' in line)  :
          line = line.lstrip()
          stripedline=line.split()
          kmxmax = float(stripedline[len(stripedline)-5])
          distance = float(stripedline[len(stripedline)-4])
          x = float(stripedline[len(stripedline)-3])
          y = float(stripedline[len(stripedline)-2])
          z = float(stripedline[len(stripedline)-1])
          logging.debug("Found Xmax altitude " + str(kmxmax) + " distance " + str(distance) + " x:" + str(x) + " y:"+ str(y) + " z:" + str(z) ) #debug level 1
          return kmxmax,distance,x,y,z
      try:
        kmxmax
      except NameError:
        logging.info('warning distance to xmax not found')
        return -1,-1,-1,-1,-1
  except:
    logging.error("GetKmXmaxFromSry:file not found or invalid:"+sry_file)
    raise
    return -1


def GetExpectedKmXmaxFromSry(sry_file,outmode="N/A"): #To do. Handle when Xmax is not found, becouse the fit didnt converge, or becouse this is not ZHAireS, or its latest version
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Exp Distance To Xmax' in line:
          line = line.lstrip()
          stripedline=line.split()
          kmxmax = float(stripedline[len(stripedline)-2])
          logging.debug("Found Expected Xmax distance " + str(kmxmax) ) #debug level 1
          return kmxmax
      try:
        kmxmax
      except NameError:
        logging.info('warning expected distance to xmax not found')
        return -1
  except:
    logging.error("GetExpectedKmXmaxFromSry:file not found or invalid:"+sry_file)
    raise
    return -1




def GetTaskNameFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Task Name:' in line:
          line = line.lstrip()
          stripedline=line.split(' ',-1)
          taskname=stripedline[len(stripedline)-1]
          taskname=taskname.replace('\n','')
          #logging.debug("Found taskname " + taskname)
          if '...' in taskname:
            base=os.path.basename(sry_file)
            taskname=os.path.splitext(base)[0]
            logging.debug("taskname contained ... (too long), using filename instead")
            print("taskname:"+taskname)

          return taskname
      try:
        taskname
      except NameError:
        logging.error('warning taskname not found, Aires has no default value, cannot continue')
        exit()
  except:
    logging.error("GetTaskNameFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetRandomSeedFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Seed of random generator:' in line:
          line = line.lstrip()
          stripedline=line.split(' ',-1)
          randomseed=stripedline[len(stripedline)-1]
          randomseed=randomseed.replace('\n','')
          return randomseed
      try:
        randomseed
      except NameError:
        logging.error('warning randomseed not found, Aires has no default value, cannot continue')
        exit()
  except:
    logging.error("GetRandomSeedFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

#output is in meters
def GetGroundAltitudeFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Ground altitude:' in line:
          line = line.lstrip()
          stripedline=line.split(':',-1)
          stripedline=stripedline[1]
          stripedline=stripedline.split(' ',-1)
          groundalt=float(stripedline[1])
          unit=stripedline[2]
          if unit == "km":
           groundalt=groundalt*1000.0
          if unit == "cm":
           groundalt=groundalt/100.0
          return groundalt
      try:
        groundalt
      except NameError:
        logging.error('warning groundalt not found, defaulting to sea level')
        return 0
  except:
    logging.error("GetGroundAltitudeFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetTimeBinFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Time Domain Bin Size:' in line:
          line = line.lstrip()
          stripedline=line.split('Time Domain Bin Size:',-1)
          stripedline=stripedline[1]
          stripedline=stripedline.lstrip()
          stripedline=stripedline.split(' ',-1)
          timebin=float(stripedline[0])
          unit=stripedline[1]
          if unit != "sec\n":
              logging.error('warning, time bin must be in seconds on the summary file, other units not supported yet')
              return -1
          return timebin*1e9

      try:
        timebin
      except NameError:
        logging.error('warning timebin not found')
        return 0
  except:
    logging.error("GetTimeBinFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetTimeWindowMinFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Antenna Time Window Min:' in line:
          line = line.lstrip()
          stripedline=line.split('Antenna Time Window Min:',-1)
          stripedline=stripedline[1]
          stripedline=stripedline.lstrip()
          stripedline=stripedline.split(' ',-1)
          timemin=float(stripedline[0])
          unit=stripedline[1]
          if unit != "sec\n":
              logging.error('warning, time window must be in seconds on the summary file, other units not supported yet')
              return -1
          return timemin*1e9

      try:
        timemin
      except NameError:
        logging.error('warning timewindow min not found')
        return 0
  except:
    logging.error("GetTimeWindowMinFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetTimeWindowMaxFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Antenna Time Window Max:' in line:
          line = line.lstrip()
          stripedline=line.split('Antenna Time Window Max:',-1)
          stripedline=stripedline[1]
          stripedline=stripedline.lstrip()
          stripedline=stripedline.split(' ',-1)
          timemax=float(stripedline[0])
          unit=stripedline[1]
          if unit != "sec\n":
              logging.error('warning, time window must be in seconds on the summary file, other units not supported yet')
              return -1
          return timemax*1e9

      try:
        timemax
      except NameError:
        logging.error('warning timewindow max not found')
        return 0
  except:
    logging.error("GetTimeWindowMaxFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetWeightFactorFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Max. stat. weight factor:' in line:
          line = line.lstrip()
          stripedline=line.split(':',-1)
          stripedline=stripedline[1]
          stripedline=stripedline.lstrip()
          stripedline=stripedline.split(' ',-1)
          wf=float(stripedline[0])
          return wf
      try:
        wf
      except NameError:
        logging.error('warning weight factor not found, default is 12')
        return 12
  except:
    logging.error("GetWeightFactorFromSry:file not found or invalid:"+sry_file)
    raise
    return -1



def GetMagneticFieldFromSry(sry_file,outmode="N/A"):
  fieldintensity=float(0.0)
  fieldinclination=float(0.0)
  fielddeclination=float(0.0)
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Geomagnetic field:' in line:
          line = line.lstrip()
          if 'Off'in line:
            fieldintensity=0.0
            fieldinclination=0.0
            fielddeclination=0.0
          else:
            stripedline=line.split('Intensity:',-1)
            intensityline=stripedline[1]
            intensityline=intensityline.lstrip()
            stripedline=intensityline.split(' ',-1)
            fieldintensity=float(stripedline[0])
        if 'I:' in line:
          line = line.lstrip()
          stripedline=line.split('I:',-1)
          inclinationline=stripedline[1]
          inclinationline=inclinationline.lstrip()
          stripedline=inclinationline.split(' ',-1)
          fieldinclination=float(stripedline[0])
          stripedline=line.split('D:',-1)
          declinationline=stripedline[1]
          declinationline=declinationline.lstrip()
          stripedline=declinationline.split(' ',-1)
          fielddeclination=float(stripedline[0])

      try:
        fieldintensity
        return fieldintensity,fieldinclination,fielddeclination
      except NameError:
        logging.error('warning MagneticField not found, defaulting to sea level')
        return 0,0,0
  except:
    logging.error("GetMagneticFieldFromSry:file not found or invalid:"+sry_file)
    raise
    return -1,-1,-1


def GetTotalCPUTimeFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Total CPU time' in line:
          line = line.lstrip()
          stripedline=line.split(':',-1)
          CPUtime=stripedline[len(stripedline)-1]
          #logging.debug("Found taskname " + taskname)
          CPUtime=CPUtime.replace('\n','')
          stripedline=CPUtime.split(' ',-1)
          time=0.0
          if(len(stripedline)>1):
            if(stripedline[-1]=="sec"):
              time=time+float(stripedline[-2])
            if(stripedline[-1]=="min"):
              time=time+60.0*int(stripedline[-2])
            if(stripedline[-1]=="hr"):
              time=time+60.0*60.0*int(stripedline[-2])
          if(len(stripedline)>3):
            if(stripedline[-3]=="sec"):
              time=time+float(stripedline[-4])
            if(stripedline[-3]=="min"):
              time=time+60.0*int(stripedline[-4])
            if(stripedline[-3]=="hr"):
              time=time+60.0*60.0*int(stripedline[-4])
          if(len(stripedline)>5):
            if(stripedline[-5]=="sec"):
              time=time + float(stripedline[-6])
            if(stripedline[-5]=="min"):
              time=time + 60.0*int(stripedline[-6])
            if(stripedline[-5]=="hr"):
              time=time + 60.0*60.0*int(stripedline[-6])
          CPUtime=time
          print("CPUTime:",CPUtime)
          return CPUtime
      try:
        CPUtime
      except NameError:
        logging.error('warning Total CPU time not found')
        return -1
  except:
    logging.error("GetTotalCPUTimeFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetHadronicModelFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Hadronic Mean Free Paths' in line:
          line = line.lstrip()
          stripedline=line.split(':',-1)
          HadronicModel=stripedline[len(stripedline)-1]
          #logging.debug("Found taskname " + taskname)
          return HadronicModel
      try:
        HadronicModel
      except NameError:
        logging.error('warning Hadronic Model not found')
        return -1
  except:
    logging.error("GetHadronicModelFromSry:file not found or invalid:"+sry_file)
    raise
    return -1


def GetRandomSeedFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Seed of random generator:' in line:
          line = line.lstrip()
          stripedline=line.split(':',-1)
          RandomSeed=stripedline[len(stripedline)-1]
          #logging.debug("Found taskname " + taskname)
          return RandomSeed
      try:
        RandomSeedl
      except NameError:
        logging.error('warning RandomSeed not found')
        return -1
  except:
    logging.error("GetRandomSeedFromSry:file not found or invalid:"+sry_file)
    raise
    return -1


def GetAiresVersionFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'This is AIRES version' in line:
          line = line.lstrip()
          stripedline=line.split('This is AIRES version',-1)
          stripedline=stripedline[-1]
          stripedline=stripedline.lstrip()
          stripedline=stripedline.split(' ',-1)
          AiresVersion=stripedline[0]
          #ogging.debug("Found Version " + AiresVersion)
          return AiresVersion
      try:
        AiresVersion
      except NameError:
        logging.error('warning Aires Version not found')
        return -1
  except:
    logging.error("GetAiresVersionFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetZHAireSVersionFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'With ZHAireS version' in line:
          line = line.lstrip()
          stripedline=line.split('With ZHAireS version',-1)
          stripedline=stripedline[-1]
          stripedline=stripedline.lstrip()
          stripedline=stripedline.split(' ',-1)
          ZHAiresVersion=stripedline[0]
          #ogging.debug("Found Version " + AiresVersion)
          return ZHAiresVersion
      try:
        ZHAiresVersion
      except NameError:
        logging.error('warning ZHAireS Version not found')
        return -1
  except:
    logging.error("GetZHAireSVersionFromSry:file not found or invalid:"+sry_file)
    raise
    return -1


def GetAtmosphericModelFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Atmospheric model' in line:
          line = line.lstrip()
          stripedline=line.split(':',-1)
          AtmosphericModel=stripedline[len(stripedline)-1]
          #logging.debug("Found taskname " + taskname)
          return AtmosphericModel
      try:
        AtmosphericModel
      except NameError:
        logging.error('warning Atmospheric Model not found')
        return -1
  except:
    logging.error("GetAtmosphericModelFromSry:file not found or invalid:"+sry_file)
    raise
    return -1


def GetSiteFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Site:' in line:
          line = line.lstrip()
          stripedline=line.split(':',-1)
          Site=stripedline[len(stripedline)-1]
          #logging.debug("Found taskname " + taskname)
          return Site
      try:
        Site
      except NameError:
        logging.error('warning Site not found')
        return -1
  except:
    logging.error("GetSiteFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetLatLongFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if '(Lat:' in line:
          line = line.lstrip()
          stripedline=line.split(':',-1)
          Lat= stripedline[1]
          Lat= Lat.lstrip()
          Lat= Lat.split(" ",-1)
          Lat= Lat[0]
          Long=stripedline[2]
          Long= Long.lstrip()
          Long= Long.split(" ",-1)
          Long= Long[0]
          return Lat,Long
      try:
        Lat
        Long
      except NameError:
        logging.error('warning Latitude or Longitude not found')
        return -1,-1
  except:
    logging.error("GetLatLongFromSry:file not found or invalid:"+sry_file)
    raise
    return -1,-1

def GetDateFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Date:' in line:
          line = line.lstrip()
          stripedline=line.split(':',-1)
          Date=stripedline[len(stripedline)-1]
          #logging.debug("Found taskname " + taskname)
          return Date
      try:
        Date
      except NameError:
        logging.error('warning Date not found')
        return -1
  except:
    logging.error("GetDateFromSry:file not found or invalid:"+sry_file)
    raise
    return -1


def GetInjectionAltitudeFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if 'Injection altitude' in line:
          line = line.lstrip()
          stripedline=line.split(':',-1) #since somtimes we can have the (D) of default and sometimes not i make sure this is whatever is after the :
          stripedline=stripedline[1]
          stripedline=stripedline.split(' ',-1)
          injalt=float(stripedline[1])
          unit=stripedline[2]
          if unit == "km":
           injalt=injalt*1000.0
          if unit == "cm":
           injalt=injalt/100.0
          return injalt
      try:
        injalt
      except NameError:
        logging.error('warning Injection Altitude not found')
        return -1
  except:
    logging.error("GetInjectionAltitudeFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def GetEnergyFractionInNeutrinosFromSry(sry_file,outmode="N/A"):
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      founditonce=0
      for line in datafile:
        if 'Neutrinos:' in line:
          if(founditonce==1):
            line = line.lstrip()
            stripedline=line.split('Neutrinos:',-1)
            stripedline=stripedline[1]
            stripedline=stripedline.lstrip()
            stripedline=stripedline.split(' ',-1)
            energyf=float(stripedline[0])
            return energyf
          founditonce=founditonce+1
      try:
        energyf
      except NameError:
        logging.error('warning NeutrinoEnergyFraction not found')
        return -1
  except:
    logging.error("GetEnergyFractionInNeutrinosFromSry:file not found or invalid:"+sry_file)
    raise
    return -1

def get_antenna_t0(xant,yant,hant, azimuthdeg, zenithdeg):
    #this code is copied from zhaires fieldinit
    #returns the t0 of the antenna in ns
    #x,y,hant is antenna position in zhaires reference frame, hant is the altitude above ground (its making flat earth asumptions for now)
    #azimuth and zenith are in ZHAireS,degrees
    cspeed = 299792458.0
    #get incoming azimut in radians
    phidirrad=azimuthdeg*np.pi/180.0
    #set phidirrad in the rangle [0,2pi]
    if(phidirrad < 0.0):
       phidirrad=phidirrad+2.0*np.pi

    zenithrad=zenithdeg*np.pi/180.0
    # Auxiliary variables
    coszenith=np.cos(zenithrad)
    sinzenith=np.sin(zenithrad)
    tanzenith=np.tan(zenithrad)


    # distance to the axis
    Rant=np.sqrt(xant*xant+yant*yant)
    #phi
    phiantrad=np.arctan2(yant,xant)
    if(phiantrad  < 0):
      phiantrad=2.0*np.pi+phiantrad

    #Adjusting time window
    #phi of antenna
    #projection of r of antenna on shower phi direction
    angle=np.absolute(phidirrad-phiantrad)

    rproj=Rant*np.cos(angle)

    dtna=(hant/coszenith + (rproj-hant*tanzenith)*sinzenith)/cspeed

    return dtna*1.0e9

def GetAntennaInfoFromSry(sry_file,outmode="N/A"):
  
  AntennaID=[]
  AntennaX=[]
  AntennaY=[]
  AntennaZ=[]
  AntennaT=[]
  Read=False
  ReadLegacy=False
  AntennaN=0
  try:
    datafile=open(sry_file,'r')
    with open(sry_file, "r") as datafile:
      for line in datafile:
        if(Read):
          stripedline=line.split()
          if(len(stripedline)==6):
            AntennaID.append(stripedline[1])
            AntennaX.append(stripedline[2])
            AntennaY.append(stripedline[3])
            AntennaZ.append(stripedline[4])
            if(stripedline[5]=="**********"):
              print("trying to recover an antenna t0")
              azimuthdeg=GetAzimuthAngleFromSry(sry_file,outmode="AIRES")
              zenithdeg=GetZenithAngleFromSry(sry_file,outmode="AIRES")
              xant=float(stripedline[2])
              yant=float(stripedline[3])
              zant=float(stripedline[4])
              ground=GetGroundAltitudeFromSry(sry_file)
              hant=zant-ground
              print(xant,yant,zant,hant,azimuthdeg,zenithdeg)
              stripedline[5]=str(get_antenna_t0(xant,yant,hant, azimuthdeg, zenithdeg))
              print(stripedline[5])


            AntennaT.append(stripedline[5])
          else:
            Read=False

            #now, i need to make the AntennaID Unique, so that i can store them in the file
            dups = {}

            for i, val in enumerate(AntennaID):
                if val not in dups:
                    # Store index of first occurrence and occurrence value
                    dups[val] = [i, 1]
                else:
                    # Special case for first occurrence
                    if dups[val][1] == 1:
                        AntennaID[dups[val][0]] += str(dups[val][1])

                    # Increment occurrence value, index value doesn't matter anymore
                    dups[val][1] += 1

                    # Use stored occurrence value
                    AntennaID[i] += str(dups[val][1])
                    
            return AntennaID,AntennaX,AntennaY,AntennaZ,AntennaT

        if(ReadLegacy):
          stripedline=line.split()
          if(len(stripedline)==5):
            AntennaX.append(stripedline[1])
            AntennaY.append(stripedline[2])
            AntennaZ.append(stripedline[3])
            AntennaT.append(stripedline[4])
            AntennaID.append("Antenna"+str(AntennaN))
            AntennaN=AntennaN+1
          else:
            ReadLegacy=False
            return AntennaID,AntennaX,AntennaY,AntennaZ,AntennaT

        elif 'Antenna|      Label      |' in line:
          Read=True
        elif 'Antenna|   X [m]' in line:
          ReadLegacy=True

  except:
    logging.error("GetAntennaInfoFromSry:file not found or invalid:"+sry_file)
    raise
    return -1,-1,-1,-1,-1

def GetLongitudinalTable(Path,TableNumber,Slant=True,Precision="Double",TaskName="Unknown"):
    #todo: check against a list of valid tables
    deletefile=False
    
    if(TaskName=="Unknown"):
      sryfile=glob.glob(Path+"/*.sry")
      idffile=glob.glob(Path+"/*.idf")
      inpfile=glob.glob(Path+"/*.inp")
      tablefile=glob.glob(Path+"/*.t"+str(TableNumber))

      if(len(tablefile)==0 and len(idffile)==0):
        logging.error("The requested table was not found, and the idf file is not present. Cannot get the table")
        return -1
    else :
      sryfile=[Path+"/"+TaskName+".sry"]
      idffile=[Path+"/"+TaskName+".idf"]
      inpfile=[Path+"/"+TaskName+".inp"]
      tablefile=glob.glob(Path+"/"+TaskName+".t"+str(TableNumber))      
          

    if(len(idffile)==1 and len(tablefile)==0):
      logging.info("could not find the table, trying to get it from the idf")
      base=os.path.basename(idffile[0])
      taskname=os.path.splitext(base)[0]

      if(len(Path)+len(taskname)<60 and len(taskname)<60):
        if(Slant==True):
            cmd=AiresPath+"/AiresExport -O a "+Path+"/"+taskname+" "+str(TableNumber)
        elif(Slant==False):
            cmd=AiresPath+"/AiresExport "+Path+"/"+taskname+" "+str(TableNumber)
        else:
            logging.error("unrecognized Slant value, please state either True/False")
            return -1
        logging.debug("abut to run "+ cmd)         
        os.system(cmd)
        tablefile=glob.glob(Path+"/*.t"+str(TableNumber))

        if(len(tablefile)==1):
            logging.debug("Table exported successfully")
            deletefile=True

      elif(len(taskname)<60):
        logging.debug("idf path+name is too long, need to copy the idf here")
        cmd="cp "+idffile[0]+" ."
        os.system(cmd)

        if(Slant==True):
            cmd=AiresPath+"/AiresExport -O a "+taskname+" "+str(TableNumber)
        elif(Slant==False):
            cmd=AiresPath+"/AiresExport "+taskname+" "+str(TableNumber)
        else:
            logging.error("unrecognized Slant value, please state either True/False")
            return -1

        os.system(cmd)
        tablefile=glob.glob("*.t"+str(TableNumber))
        logging.debug(tablefile)

        if(len(tablefile)==1):#this means we are not wo
            logging.debug("Table exported successfully")
            deletefile=True
            cmd="rm "+taskname+".idf"
            os.system(cmd)
            cmd="rm "+taskname+".lgf"
            os.system(cmd)
      else:
       logging.debug("task name is to long, AIRES does not support that!")

    if(len(tablefile)==1):
      logging.debug("reading file")

      from numpy import loadtxt

      if(Precision=="Double"):
        numpyarray=loadtxt(tablefile[0],usecols=(1,2),dtype='f8')
      elif(Precision=="Simple"):
        numpyarray=loadtxt(tablefile[0],usecols=(1,2),dtype='f4')
      else:
        logging.error("unrecognized precison, please state either Double or Simple")
        return -1

      if(deletefile==True):
        cmd="rm "+tablefile[0]
        os.system(cmd)
        logging.debug("Table deleted successfully")

      return numpyarray
    else:
      logging.error("The requested table was not found and could not be regenerated. Sorry")
      return -1


def GetLateralTable(Path,TableNumber,Density=True,Precision="Double"):
    #todo: check against a list of valid tables
    deletefile=False
    sryfile=glob.glob(Path+"/*.sry")
    idffile=glob.glob(Path+"/*.idf")
    inpfile=glob.glob(Path+"/*.inp")
    tablefile=glob.glob(Path+"/*.t"+str(TableNumber))

    if(len(tablefile)==0 and len(idffile)==0):
      logging.error("LDF.The requested table was not found, and the idf file is not present. Cannot get the table")
      return -1

    if(len(idffile)==1 and len(tablefile)==0):
      logging.info("LDF.could not find the table, trying to guess it from the idf")
      base=os.path.basename(idffile[0])
      taskname=os.path.splitext(base)[0]

      if(len(Path)+len(taskname)<60 and len(taskname)<60):

          if(Density==True):
            cmd=AiresPath+"/AiresExport -O dX "+Path+"/"+taskname+" "+str(TableNumber)
          elif(Density==False):
            cmd=AiresPath+"/AiresExport -O X "+Path+"/"+taskname+" "+str(TableNumber)
          else:
            logging.error("LDF.unrecognized Density value, please state either True/False")
            return -1

          os.system(cmd)
          tablefile=glob.glob(Path+"/*.t"+str(TableNumber))

          if(len(tablefile)==1):
            logging.debug("LDF.Table exported successfully")
            deletefile=True

      elif(len(taskname)<60):
          logging.debug("LDF. idf path+name is too long, need to copy the idf here")
          cmd="cp "+idffile[0]+" ."
          os.system(cmd)

          if(Density==True):
            cmd=AiresPath+"/AiresExport -O dX "+taskname+" "+str(TableNumber)
          elif(Density==False):
            cmd=AiresPath+"/AiresExport -O X "+taskname+" "+str(TableNumber)
          else:
            logging.error("LFD.unrecognized Density value, please state either True/False")
            return -1

          os.system(cmd)
          tablefile=glob.glob("*.t"+str(TableNumber))

          if(len(tablefile)==1):
            logging.debug("LDF.Table exported successfully")
            deletefile=True
            cmd="rm "+taskname+".idf"
            os.system(cmd)
            cmd="rm "+taskname+".lgf"
            os.system(cmd)

      else:
         logging.debug("LDF.task name is to long, AIRES does not support that!")

    if(len(tablefile)==1):
      logging.debug("LDF.reading file")

      from numpy import loadtxt

      if(Precision=="Double"):
        numpyarray=loadtxt(tablefile[0],usecols=(1,2),dtype='f8')
      elif(Precision=="Simple"):
        numpyarray=loadtxt(tablefile[0],usecols=(1,2),dtype='f4')
      else:
        logging.error("LDF.unrecognized precison, please state either Double or Simple")
        return -1

      if(deletefile==True):
        cmd="rm "+tablefile[0]
        os.system(cmd)
        logging.debug("LDF.Table deleted successfully")

      return numpyarray
    else:
      logging.error("LDF.The requested table was not found and could not be regenerated. Sorry")
      return -1






#this gets the effective refraction index from poitn R0 to xant,yant,zant (default to the core position), all in meters, but kr in 1/km
#ns is the REFRACIVITY AR SEA LEVEL, kr is the scale height in km.
import numpy as np
#NOTE: Use the GetZHSEffectiveactionIndex below to get exactly the same results than in zhaires
# ns= refractivity at sea levele
# kr= scale height

def GetEffectiveRefractionIndex(x0,y0,z0,ns=325,kr=-0.128,zant=0,xant=0,yant=0,stepsize = 20000): #NOTE THAT THE ORDER OF ANTENNA POSITIONS IS ZANT,XANT,YANT (historic compatibility reasons)
        rearth=6370949.0
        R02=x0*x0+y0*y0  #notar que se usa R02, se puede ahorrar el producto y la raiz cuadrada
        h0=(np.sqrt((z0+rearth)*(z0+rearth) + R02 ) - rearth)/1.E3 #altitude of emission, in km

        rh0 = ns*np.exp(kr*h0) #refractivity at emission
        n_h0=1.E0+1.E-6*rh0 #n at emission
#        print("n_h0",n_h0,ns,kr,x0,y0,z0,h0,rh0)

        hd=(zant)/1.E3 #detector altitude

#       Vector from detector to average point on track. Making the integral in this way better guaranties the continuity
#       since the choping of the path will be always the same as you go farther away. If you start at your starting point, for a given geometry,
#       the choping points change with each starting position.

        ux = x0-xant
        uy = y0-yant         #the antenna position, considered to be at the core
        uz = z0-zant

        Rd=np.sqrt(ux*ux + uy*uy)
        kx=ux/Rd
        ky=uy/Rd #!k is a vector from the antenna to the track, that when multiplied by Rd will end in the track and sumed to antenna position will be equal to the track positon
        kz=uz/Rd

#       integral starts at ground
        nint=0
        sum=0.E0

        currpx=0+xant
        currpy=0+yant    #!current point (1st antenna position)
        currpz=zant
        currh=hd
        
        print(Rd)

        while(Rd > stepsize): #if distance projected on the xy plane is more than 10km
          nint=nint+1
          nextpx=currpx+kx*stepsize
          nextpy=currpy+ky*stepsize           #this is the "next" point
          nextpz=currpz+kz*stepsize

          nextR2=nextpx*nextpx + nextpy*nextpy #!se usa el cuadrado, se puede ahorrar la raiz cuadrada
          nexth=(np.sqrt((nextpz+rearth)*(nextpz+rearth) + nextR2) - rearth)/1.E3

          if(np.absolute(nexth-currh) > 1.E-10  ):   #check that we are not going at constant height, if so, the refraction index is constant
              sum=sum+(np.exp(kr*nexth)-np.exp(kr*currh))/(kr*(nexth-currh))
          else:
              sum=sum+np.exp(kr*currh)

          currpx=nextpx
          currpy=nextpy
          currpz=nextpz  #Set new "current" point
          currh=nexth

          Rd=Rd-stepsize #reduce the remaining lenght
        #enddo

        #when we arrive here, we know that we are left with the last part of the integral, the one closer to the track (and maybe the only one)

        nexth=h0

        if(np.absolute(nexth-currh) > 1.E-10 ): #check that we are not going at constant height, if so, the refraction index is constant
          sum=sum+(np.exp(kr*nexth)-np.exp(kr*currh))/(kr*(nexth-currh))
        else:
          sum=sum+np.exp(kr*currh)

        nint=nint+1
        avn=ns*sum/nint
        n_eff=1.E0+1.E-6*avn #average (effective) n
        return n_eff

#This is exactly what is computed in ZHAires
def GetZHSEffectiveRefractionIndex(x0,y0,z0,xant=0,yant=0,zant=0,ns=325,kr=-0.1218,stepsize = 20000):
#      rearth=6371007.0 #new aires
      rearth=6370949.0 #19.4.0
#     Variable n integral calculation ///////////////////
      R02=x0*x0+y0*y0  #!notar que se usa R02, se puede ahorrar el producto y la raiz cuadrada (entro con injz-zXmax -> injz-z0=zXmax
      h0=(np.sqrt( (z0+rearth)*(z0+rearth) + R02 ) - rearth)/1E3    #!altitude of emission

      rh0=ns*np.exp(kr*h0) #!refractivity at emission (this
      n_h0=1+1E-6*rh0 #!n at emission
#        write(*,*) "n_h0",n_h0,ns,kr,x0,y0,injz-z0,h0,rh0
      modr=np.sqrt(R02)

#      print("ZHS n_h0",n_h0,ns,kr,x0,y0,z0,h0,rh0)
        
      if(modr > 1000): #! if inclined shower and point more than 20km from core. Using the core as reference distance is dangerous, its invalid in upgoing showers

#         Vector from average point of track to observer.
          ux = xant-x0
          uy = yant-y0
          uz = zant-z0

#         divided in nint pieces shorter than 10km
          nint=int((modr/stepsize)+1)
          kx=ux/nint
          ky=uy/nint       #k is vector from one point to the next
          kz=uz/nint
#
          currpx=x0
          currpy=y0        #current point (1st is emission point)
          currpz=z0
          currh=h0
#
          sum=0
          for iii in range(0,nint):
            nextpx=currpx+kx
            nextpy=currpy+ky #!this is the "next" point
            nextpz=currpz+kz
            nextR2=nextpx*nextpx + nextpy*nextpy
            nexth=(np.sqrt((nextpz+rearth)*(nextpz+rearth) + nextR2) - rearth)/1E3
#c
            if(np.abs(nexth-currh) > 1E-10  ):
              sum=sum+(np.exp(kr*nexth)-np.exp(kr*currh))/(kr*(nexth-currh))
            else:
              sum=sum+np.exp(kr*currh)
#            endif
#c
            currpx=nextpx
            currpy=nextpy
            currpz=nextpz  #!Set new "current" point
            currh=nexth
#c
          avn=ns*sum/nint
#          print*,"avn:",avn
          n_eff=1+1E-6*avn #!average (effective) n
#c
      else:
#c         withouth integral
          hd=zant/1E3 #!detector altitude(flat earth!)
#
          if(np.abs(hd-h0) > 1E-10):
            avn=(ns/(kr*(hd-h0)))*(np.exp(kr*hd)-np.exp(kr*h0))
#            print*,"avn2:",avn
          else:
            avn=ns*np.exp(kr*h0)
#           print*,"avn3:",avn
#            print *,"Effective n: h0=hd"
#          endif
          n_eff=1+1E-6*avn #!average (effective) n
#        endif
#c     ///////////////////////////////////////////////////
      return n_eff
#      end



def ReadAiresSry(sry_file,outmode="N/A"):

  zen=GetZenithAngleFromSry(sry_file,outmode)
  azim=GetAzimuthAngleFromSry(sry_file,outmode)
  energy=GetEnergyFromSry(sry_file,outmode)
  primary=GetPrimaryFromSry(sry_file,outmode)
  xmax=GetSlantXmaxFromSry(sry_file,outmode)
  kmxmax,distance,x,y,z=GetKmXmaxFromSry(sry_file,outmode)
  taskname=GetTaskNameFromSry(sry_file,outmode)
  return zen,azim,energy,primary,xmax,distance,taskname

def ReadAiresLgf(lgf_file,outmode="N/A"):

  zen=GetZenithAngleFromSry(lgf_file,outmode)
  azim=GetAzimuthAngleFromSry(lgf_file,outmode)
  energy=GetEnergyFromSry(lgf_file,outmode)
  primary=GetPrimaryFromSry(lgf_file,outmode)
  taskname=GetTaskNameFromSry(lgf_file,outmode)
  return zen,azim,energy,primary,-1,-1,taskname


def GetStatusFromStatus(status_file):
  try:
    datafile=open(status_file,'r')
    with open(status_file, "r") as datafile:
      for line in datafile:
        if 'Aires_Msg' in line:
          line = line.lstrip()
          stripedline=line.split('=',-1)
          status=stripedline[len(stripedline)-1]
          status=status.replace('\'','')
          status=status.replace('\n','')
          #logging.debug("Found status " + status)
          return status
      try:
        status
      except NameError:
        logging.error('warning status not found in status file')
        return 'not found'
  except:
    logging.error("GetStatusFromStatus:file not found or invalid:"+status_file)
    return -1

def GetTmpFromDirs(dirs_file):
  try:
    datafile=open(dirs_file,'r')
    with open(dirs_file, "r") as datafile:
      for line in datafile:
        if 'Aires_DRandomfn' in line:
          line = line.lstrip()
          stripedline=line.split('=',-1)
          tmp=stripedline[len(stripedline)-1]
          tmp=tmp.replace('\'','')
          tmp=tmp.replace('\n','')
          #logging.debug("Found Tmp " + tmp)
          return tmp
      try:
        tmp
      except NameError:
        logging.error('warning Aires_DRandomfn not found in dirs file')
        return 'not found'
  except:
    logging.error("GetTmpFromSDirs:file not found or invalid:"+dirs_file)
    return -1

def DeprecatedReadAiresSry(sry_file,outmode="GRAND"):

    try:
     datafile=open(sry_file,'r')

    except:
      logging.error("ReasAiresSry:file not found or invalid:"+sry_file)
      raise
      return -1, -1, -1, -1, -1

    for line in datafile:
        ## print(line) #debug level 2

        if 'Primary zenith angle:' in line:
            line = line.lstrip()
            stripedline=line.split(' ',-1)
            zen=float(stripedline[len(stripedline)-2])
            if outmode == 'GRAND':
                zen = 180-zen  #conversion to GRAND convention i.e. pointing towards antenna/propagtion direction
            logging.debug('Found Zenith ' + str(zen))

        if 'Primary azimuth angle:' in line:
            line = line.lstrip()
            stripedline=line.split(' ',-1)
            azim = float(stripedline[len(stripedline)-2])

            if outmode == 'GRAND':
                azim=azim+180 #conversion to GRAND convention i.e. pointing towards antenna/propagtion direction
                if azim>=360:
                    azim= azim-360
            logging.debug('Found Azimuth ' + str(azim))

        if 'Primary energy:' in line:
            line = line.lstrip()
            stripedline=line.split(' ',-1)
            try:
              energy
            except NameError:
              energy = float(stripedline[len(stripedline)-2])
              unit= str(stripedline[len(stripedline)-1])

              if outmode == 'GRAND':
                if unit == "eV\n":
                    energy = energy *1e-18
                if unit == "KeV\n":
                    energy = energy *1e-15
                if unit == "MeV\n":
                    energy = energy *1e-12
                if unit == "GeV\n":
                    energy = energy *1e-9
                if unit == "TeV\n":
                    energy = energy *1e-6
                if unit == "PeV\n":
                    energy = energy *1e-3
                if unit == "EeV\n":
                    energy = energy

              if outmode == 'AIRES':
                if unit == "eV\n":
                    energy = energy *1e-9
                if unit == "KeV\n":
                    energy = energy *1e-6
                if unit == "MeV\n":
                    energy = energy *1e-3
                if unit == "GeV\n":
                    energy = energy
                if unit == "TeV\n":
                    energy = energy *1e3
                if unit == "PeV\n":
                    energy = energy *1e6
                if unit == "EeV\n":
                    energy = energy *1e9

              logging.debug('Found Energy ' + str(energy)) #debug level 1

        if 'Primary particle:' in line:
            line = line.lstrip()
            stripedline=line.split(' ',-1)
            if(len(stripedline)==3):
              primarytype = str(stripedline[len(stripedline)-1])
            elif(len(stripedline)==5):
              primarytype = str(stripedline[len(stripedline)-3])
            elif(len(stripedline)==6):
              primarytype = str(stripedline[len(stripedline)-4])
            elif(len(stripedline)==7):
              primarytype = str(stripedline[len(stripedline)-5])
            else:
              primarytype = "unknown"

            logging.debug('Found Primary ' + primarytype) #debug level 1


        if 'Sl. depth of max. (g/cm2)' in line:
            line = line.lstrip()
            stripedline=line.split(' ',-1)
            xmax = float(stripedline[len(stripedline)-1])
            logging.debug('Found Xmax ' + str(xmax)) #debug level 1

    try:
        zen
    except NameError:
        zen = 0 #If no zenith angle was included in the input file, AIRES defaults to 0
        if outmode == 'GRAND':
          zen = 180-0 #that translates to 180 in GRAND
    try:
        azim
    except NameError:
        azim = 0 #If no azimuth angle was included in the input file, AIRES defaults to 0
        if outmode == 'GRAND':
          azim = 0+180 # that translates to 180
    try:
        energy
    except NameError:
        logging.error('warning energy not found, Aires has no default value,  cannot continue')
        exit()
    try:
        primarytype
    except NameError:
        logging.error('warning primary not found, Aires has no default value, cannot continue')
        exit()

    try:
        xmax
    except NameError:
        logging.info('warning xmax not found')
        xmax=-1

    return zen,azim,energy,primarytype,xmax



if __name__ == '__main__':

    path = sys.argv[1]
    outmode = 'AIRES'
    print(ReadAiresSry(path,outmode))