#======================================================================================
# This script is an example showing how to use the LineAnalysisScripting module with :
# - two molecules (12CO & 13CO)                                                         
# - chi2 computation in a regular grid                                               
# - one component with the LTE model                                                  
#
# Data file in myDirInput : I16293_CO.fus
# Files created in myDirOutput : myName_bestModel.lis, myName.lam, myName.dat, myName.png
# (see below for the definitions of myDirInput, myDirOutput and myName)
# 
# To run this script as it is, the only thing you need to specify is the path to your python executable (myPython variable).
#
# Last update: 12-jan-2026
# Tested on: CASSIS 6.5
#
# CASSIS scripts and scripting documentation :
#   - http://cassis.irap.omp.eu/docs/CASSIS_scripts_summary.pdf
#   - http://cassis.irap.omp.eu/docs/script/README.html
#======================================================================================

import os
import time
import ScriptEnvironment
import subprocess
from Range import Range
from LineAnalysisScripting import UserInputs
from Component import Component
from eu.omp.irap.cassis.properties import Software
from java.io import File



# INPUTS

# Define the frequency ranges for the lines
Range.unit = "GHz"  # Possible units are GHz and km/s, MHz, cm-1 and micrometer
r1 = Range(550.91, 550.94)  # or v1 = Range(-5, 10.0) if Range.unit = "km/s"
r2 = Range(576.25, 576.29)  # or v2 = Range(-5, 10.0) if Range.unit = "km/s"
r3 = Range(1101.32, 1101.37)  # or v2 = Range(-5, 10.0) if Range.unit = "km/s"
r4 = Range(1151.95, 1152.0)  # or v2 = Range(-5, 10.0) if Range.unit = "km/s"

# Set some useful variables
sourceName = "I16293"
speciesName = "12CO_13CO"
myModel = "1comp_lte"
myName =  sourceName + "_" + speciesName + "_" + myModel

# Change the lines below for your own myDirInput & myDirOutput folder for your own scripts
myDirInput = Software.getCassisPath() + "/delivery/script/data/"
myDirOutput = Software.getCassisPath() + "/delivery/script/outputs/"
if not os.path.exists(myDirOutput):
    os.makedirs(myDirOutput)

inputFile = myDirInput + "I16293_CO.fus"
outputFile = myDirOutput + myName + ".dat"

# For contour plots of chi2 values, set the path for your python executable below 
# (see requirements.txt in CASSIS_PATH/delivery/script/python/triangle_plots/ for the required packages)
myPython = "~/PycharmProjects/cassis-scripts/.venv/bin/python3"

myPythonScript = Software.getCassisPath() + "/delivery/script/python/triangle_plots/plot_chi2_RG.py"
# For more information on plot_chi2_RG.py, please open the file and read the explanation at the top of the script. 

# CONFIG

userInputs  = UserInputs(
    inputFile = inputFile,
    outputFile  = outputFile,
    #telescope = {apex: [0, 1], hifi: [2, 3]},
    telescope = "hifi",
    tuningRange = [550.0, 1152.0],  # in GHz
    tuningBand  = 60.0,  # in km/s
    aijMin = 0.0,
    eup = [0.0, 350.0],
    kup = ["*","*"],
    template = "ISM",  # or "Full CDMS" or "Full JPL" or "Full VASTEL" etc...
    moltags = [28503, 29501],
    tmb2ta = True,
    isoUnique = False,
    plotTitle = myName,
    warning = True,
    continuum = "continuum-0",  # file containing the continuum adapted to your observations
    # continuum = 3.4,  # or value of the continuum adapted to your observations
    observing_mode = "PSw/DBSw",  #or observing_mode = "FSw",

    # Enter here the lines and the corresponding ranges to be taken into account in the computation.
    # The lines are sorted by frequency.
    selectedLines = {"1": r1, "2": r2, "3": r3, "4": r4},
    # if selection on a velocity range, edit the line below
    # selectedLines = {"1": v1, "2": v2, "3": v3, "4": v4},

    # rmsLines data (in K)
    rmsLines = {"1-4": 0.200},

    # calibration for the lines
    calLines = {"1-4": 0.0}
)

# MODEL INPUTS

# Parameters for the first component
comp_1 = Component(
    model = "lte", #  or: model = "radex",
    # Parameters needed for LTE or RADEX model
    nmol = {'min': 9.5e15, 'max': 1.05e16, 'nstep': 10, 'log_mode': True},
    temp = {'min': 95.0, 'max':105.0, 'nstep':10, 'log_mode': False},
    fwhm = {'min': 5.0, 'max':5.0, 'nstep':6, 'log_mode': False},
    size = {'min': 33333.0, 'max':33333.0, 'nstep': 3, 'log_mode': False},
    vlsr = {'min': 3.80, 'max':3.80, 'nstep': 3, 'log_mode': False},
    iso = {'min': 40, 'max':60, 'nstep': 10, 'log_mode': False},  # if moltags are two isotopologues
    interacting = True,
    # the following next few lines are only needed for RADEX
    collisioner = ["p_H2", "o_H2"],
    n_p_H2 = {'min': 1.0e4, 'max': 1.0e8, 'nstep': 2, 'log_mode': True},  # density for the collisioner(s)
    n_o_H2 = {'min': 1.0e4, 'max': 1.0e8, 'nstep': 2, 'log_mode': True},  # density for the collisioner(s)
    collisionFile = ["co.dat", "13co.dat"],  # collisional files in the CASSIS/database/lamda directory
    geometry = "sphere" # or: geometry= "slab"
)

# Parameters for the second component
comp_2 = Component(
    model="lte",  # or: model = "radex",
    # Needed for LTE or RADEX model
    nmol = {'min': 1.0e17, 'max': 3.0e17, 'nstep': 2, 'log_mode': True},
    temp = {'min': 100.0, 'max': 200.0, 'nstep': 2, 'log_mode': False},
    fwhm = {'min': 3.5, 'max': 3.5, 'nstep': 4, 'log_mode': False},
    size = {'min': 30.0, 'max': 30.0, 'nstep': 3, 'log_mode': False},
    vlsr = {'min': 3.80, 'max': 3.80, 'nstep': 3, 'log_mode': False},
    iso = {'min': 10, 'max': 100, 'nstep': 3, 'log_mode': False},  # if moltags are two isotopologues
    interacting = True,
    # the following next few lines are only needed for RADEX
    collisioner = ["p_H2", "o_H2"],
    n_p_H2 = {'min': 1.0e4, 'max': 1.0e8, 'nstep': 2, 'log_mode': True},  # density for the collisioner(s)
    n_o_H2 = {'min': 1.0e4, 'max': 1.0e8, 'nstep': 2, 'log_mode': True},  # density for the collisioner(s)
    collisionFile = ["co.dat", "13co.dat"],  # collisional files in the CASSIS/database/lamda directory
    geometry = "sphere" # or: geometry= "slab"
)

# Execution time beginning
timeStart = time.time()

# Compute the chi2 min using the regular grid and write the data
userInputs.computeChi2MinUsingRG(comp_1)  # if one component
# userInputs.computeChi2MinUsingRG(comp_1,comp_2)  # if two components

# CASSIS Execution time ending
timeEnd = time.time()

print("CASSIS execution time = %.3f seconds." % (timeEnd - timeStart))

# SAVE AND DISPLAY THE RESULTS

# A. Plot the best model and save the corresponding spectra and config files
bestLineModel = userInputs.plotBestModel(moltag=[28503, 29501], overSampling=5)
bestLineModel.saveConfig(File(myDirOutput + myName+".lam"))
bestPhysicalModels  = userInputs.getBestPhysicalModels()
userInputs.saveBestPhysicalModels(myDirOutput + myName + "_bestModel.lis")

# B. Contour plots for chi2
trianglePlot = [myPython + " " + myPythonScript + " " + outputFile]
subprocess.Popen(trianglePlot, shell=True)
