#======================================================================================
# This script is an example showing how to use the LineAnalysisScripting module with : 
# - one molecule (o-H2CO)
# - chi2 computation with a regular grid
# - one components with the radex model
#
# Data files in myDirInput : Orion_Synthetic_H2CO.lis, Orion_Synthetic_Continuum
# 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 Component import Component
from LineAnalysisScripting import UserInputs
from eu.omp.irap.cassis.properties import Software
from java.io import File


# INPUTS

# Define the velocity ranges for the lines
Range.unit = "km/s"  # Possible units are GHz and km/s, MHz, cm-1 and micrometer
v1 = Range(0.0, 12.0)

# Set some useful variables
sourceName = "Orion_Synthetic"
speciesName = "o-H2CO"
myModel = "2comp_radex"
myName = sourceName + "_" + speciesName + "_" + myModel

# Change the two 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+"Orion_Synthetic_H2CO.lis"
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,
    telescope = "hifi",
    tuningRange = [480, 640],
    tuningBand = 60,
    aijMin = 1.0e-4,
    eup = [0.0, 250.0],
    template = "All Species",
    moltags = [30591],  # ortho-H2CO
    tmb2ta = False,
    isoUnique = False,
    outputFile = outputFile,
    plotTitle = myName,
    warning = True,

    # Enter here the lines and the corresponding ranges to be taken into account in the computation
    # The lines should be sorted by increasing frequency.
    selectedLines = {"1,4-9": [v1]},

    # rms of the data around each selected line (in K)
    rmsLines = {"1,4-9": [0.025]},

    # calibration accuracy of the data around each selected line
    calLines = {"1,4-9": [0.15]},

    continuum = myDirInput + "Orion_Synthetic_Continuum",
)

# MODEL INPUTS

comp_1 = Component(
    model = "radex",
    nmol = {'min': 1.1e14, 'max': 1.1e14, 'nstep': 5, 'log_mode': True},
    temp = {'min': 90.0, 'max': 90.0, 'nstep': 10, 'log_mode': False},
    # warning: for radex, use a T range within the T range valid for the collision file
    fwhm = {'min': 5.0, 'max': 5.0, 'nstep': 1, 'log_mode': False},
    size = {'min': 10, 'max': 10, 'nstep': 1, 'log_mode': False},
    vlsr = {'min': 6.2, 'max': 6.2, 'nstep': 1, 'log_mode': False},
    iso = {'min': 70, 'max': 70, 'nstep': 1, 'log_mode': False},  # if moltags are two isotopologues
    interacting = True,
    # Parameters below are only needed for RADEX
    collisioner = ["o_H2"],
    n_o_H2 = {'min': 1.0e8,  'max': 1.0e8,  'nstep': 5,  'log_mode': True},  # density for the collisioner(s)
    collisionFile = ["o-H2CO-o-p-H2-cdms.dat"],
    geometry = "sphere"
)

comp_2 = Component(
    model = "radex",
    nmol = {'min': 1.0e14, 'max': 5.0e14, 'nstep': 7,  'log_mode': True},
    temp = {'min': 22.0, 'max': 32.0, 'nstep': 7, 'log_mode': False},
    # warning: for radex, use a T range within the T range valid for the collision file
    fwhm = {'min': 1.0, 'max': 1.0, 'nstep': 1, 'log_mode': False},
    size = {'min': 10, 'max': 10, 'nstep': 1, 'log_mode': False},
    vlsr = {'min': 7.0, 'max': 7.0, 'nstep': 1, 'log_mode': False},
    iso = {'min': 70, 'max': 70, 'nstep': 1, 'log_mode': False},  # if moltags are two isotopologues
    interacting = True,
    # Parameters below are only needed for RADEX
    collisioner = ["o_H2"],
    n_o_H2 = {'min': 3.0e6, 'max': 3.0e7, 'nstep': 7, 'log_mode': True},  # density for the collisioner(s)
    collisionFile = ["o-H2CO-o-p-H2-cdms.dat"],
    geometry = "sphere"
)

# Execution time beginning
timeStart = time.time()

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

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

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

# ANALYSIS OF THE RESULTS

# A. Plot the best model and save the corresponding spectra and config files
lineModel = userInputs.plotBestModel(moltag=30591, overSampling=3, tuningBand=60, telescope="hifi", lines="1,4-9")
lineModel.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)
