#======================================================================================
# This script is an example showing how to use the LineAnalysisScripting module with:
# - one molecule (H2CO)
# - chi2 computation in a regular grid
# - continuum from file, two components with the LTE 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 ScriptEnvironment
import time
import subprocess
from Component import Component
from Range import Range
from LineAnalysisScripting import UserInputs
from eu.omp.irap.cassis.properties import Software
from java.io import File


# INPUTS

Range.unit = "km/s"  # Possible units are GHz and km/s, MHz, cm-1 and micrometer
v1 = Range(0.0, 15.0)

sourceName = "Orion_Synthetic"
speciesName = "p-H2CO"
myModel = "2comp_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)

myCont = myDirInput+"Orion_Synthetic_Continuum"
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],  # in GHz
    tuningBand = 60, # in km/s
    aijMin = 1.0e-4,
    eup = [0.0, 250.0],
    template = "Full VASTEL",  # or "CDMS" or "JPL" or "ISM" etc...
    moltags = [30581],  # para-H2CO
    tmb2ta = False,
    outputFile = outputFile,
    plotTitle = myName,
    warning = True,

    # Enter here the lines and the corresponding ranges to be selected for the computation
    selectedLines = {"1-6": [v1]},
    # Lines data
    rmsLines = {"1-6": [0.025]},
    calLines = {"1-6": [0.15]},
    continuum = myCont,
    observing_mode = "PSw/DBSw"  # or observing_mode = "FSw"
)

# MODEL INPUTS

# Parameters for the first component
comp_1 = Component(
    model = "lte",
    nmol = {'min': 3.0e13, 'max': 4.0e13, 'nstep': 5, 'log_mode': True},
    temp = {'min': 60.0, 'max': 110.0, 'nstep': 7, 'log_mode': False},
    fwhm = {'min': 5.0, 'max': 5.0, 'nstep': 4, 'log_mode': False},
    size = {'min': 10, 'max': 10, 'nstep': 3, 'log_mode': False},
    vlsr = {'min': 6.0, 'max': 6.0, 'nstep': 3, 'log_mode': False},
    iso = {'min': 100, 'max': 100, 'nstep': 3, 'log_mode': False},  # if moltags are two isotopologues
    interacting = True,
    # Parameters below are only needed for RADEX
    collisioner = ["p_H2", "o_H2"],
    n_p_H2 = {'min': 1.0e8, 'max': 1.0e8, 'nstep': 1, 'log_mode': True},  # density for the collisioner(s)
    n_o_H2 = {'min': 1.0e8, 'max': 1.0e8, 'nstep': 1, 'log_mode': True},  # density for the collisioner(s)
    collisionFile = ["p-H2CO-o-p-H2-cdms.dat"],  # collisional files in the CASSIS/database/lamda directory
    geometry = "sphere"  # or "slab"
)

# Parameters for the second component
comp_2 = Component(
    model = "lte",
    nmol = {'min': 1.0e12, 'max': 6.0e12, 'nstep': 7, 'log_mode': True},
    temp = {'min': 19.0, 'max': 28.0, 'nstep': 7, 'log_mode': False},
    fwhm = {'min': 1.0, 'max': 1.0, 'nstep': 4, 'log_mode': False},
    size = {'min': 3333, 'max': 3333, 'nstep': 4, 'log_mode': False},
    vlsr = {'min': 7.0, 'max': 7.0, 'nstep': 5, 'log_mode': False},
    iso = {'min': 100, 'max': 100, 'nstep': 3, 'log_mode': False},  # if moltags are two isotopologues
    interacting = True,
    # Parameters below are only needed for RADEX
    collisioner = ["p_H2","o_H2"],
    n_p_H2 = {'min': 1.0e8, 'max': 1.0e8, 'nstep': 1, 'log_mode': True},  # density for the collisioner(s)
    n_o_H2 = {'min': 1.0e8, 'max': 1.0e8, 'nstep': 1, 'log_mode': True},  # density for the collisioner(s)
    collisionFile = ["p-H2CO-o-p-H2-cdms.dat"],  # collisional files in the CASSIS/database/lamda directory
    geometry = "sphere"  # or "slab"
)

# Execution time beginning
timeStart = time.time()

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

# 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=30581, overSampling=3, tuningBand=60)
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)
