#!/bin/sh  
# \
exec oagtclsh "$0" "$@"

# Compute elastic gas scattering lifetime from DA statistics

if {![info exists env(OAG_TOP_DIR)]} { set env(OAG_TOP_DIR) /usr/local }
set auto_path [linsert $auto_path 0  $env(OAG_TOP_DIR)/oag/apps/lib/$env(HOST_ARCH)]

APSStandardSetup

set allowedNameList [list Min Max Median Clipped]
for {set p 5} {$p<=100} {incr p 5} {
    lappend allowedNameList ${p}Percentile
}

set usage {usage: elasticScatteringLifetime -aperture <daFile> -name <string>(10Percentile) -twiss <twissFile> -output <filename> [-sPeriodic <value>27.61] [-pressure nTorr(2)]}
set aperture ""
set twiss ""
set name 10Percentile
set output ""
set sPeriodic 27.61
set pressure 2
set args $argv
if {[APSStrictParseArguments {aperture name twiss output sPeriodic pressure}] || ![string length $aperture] || ![string length $twiss] || ![string length $output]} {
    return -code error "$usage"
}

if ![file exists $aperture] {
   return -code error "not found: $aperture"
}
if ![file exists $twiss] {
   return -code error "not found: $twiss"
}
if [file exists $output] {
   return -code error "in use: $output"
}
if [lsearch -exact $allowedNameList $name]==-1 {
    return -code error "not valid: $name"
}

set tmpRoot [APSTmpString]
APSAddToTempFileList $tmpRoot.1 $tmpRoot.2

# Find beta functions at ID straight, use only first sector's worth of data
# Sort into increasing s order, remove duplicates
# 
exec sddsprocess $twiss -pipe=out -process=beta*,first,%s0 -filter=col,s,0,$sPeriodic \
    | sddssort -pipe -column=s -unique \
    | sddsconvert -pipe -retain=col,betax,betay,s -retain=parameter,pCentral,beta?0 \
    | sddsexpand -pipe=in $tmpRoot.1
set nTw [exec sdds2stream -npage=bare $tmpRoot.1]

set nAp   [exec sdds2stream -npage=bare $aperture]
if $nAp>1 {
    eval exec sddscombine [APSReplicateItem -item $tmpRoot.1 -number $nAp] $tmpRoot.2
    set twissFile $tmpRoot.2
} else {
    set twissFile $tmpRoot.1
}


# Gas fractional pressures
set fGas(H2) 0.65
set fGas(H2O) 0.02
set fGas(CH4) 0.07
set fGas(CO) 0.22
set fGas(CO2) 0.04
set fGas(N2) 0.0

# Gas constituents (nAtomTypes=# of types of atoms, Z=nuclear charge, atomCount=# atoms of each type)
set nAtomTypes(H2) 1
set Z(H2,1) 1
set atomCount(H2,1) 2

set nAtomTypes(H2O) 2
set Z(H2O,1) 1
set atomCount(H2O,1) 2
set Z(H2O,2) 6
set atomCount(H2O,2) 1

set nAtomTypes(CH4) 2
set Z(CH4,1) 8
set atomCount(CH4,1) 1
set Z(CH4,2) 1
set atomCount(CH4,2) 4

set nAtomTypes(CO) 2
set Z(CO,1) 8
set atomCount(CO,1) 1
set Z(CO,2) 6
set atomCount(CO,2) 1

set nAtomTypes(CO2) 2
set Z(CO2,1) 8
set atomCount(CO2,1) 1
set Z(CO2,2) 6
set atomCount(CO2,2) 2

set nAtomTypes(N2) 1
set Z(N2,1) 7
set atomCount(N2,1) 2

set fSZ2Sum 0.0
foreach gas [array names fGas] {
    set f $fGas($gas)
    set sum 0.0
    for {set con 1} { $con <= $nAtomTypes($gas) } {incr con} {
        set sum [expr $sum+$Z($gas,$con)*$Z($gas,$con)*$atomCount($gas,$con)]
    }
    set fSZ2Sum [expr $fSZ2Sum+$f*$sum]
}
# Convert to Pascals
set P [expr 133.3224*1e-9*$pressure]

# The constant c1 is 2x the value in AOP-TN-2014-009 to account for the fact that we
# are only integrating phi:[0, pi], not [0, 2*pi]

eval exec sddscombine [APSReplicateItem -number $nTw -item $aperture] -pipe=out -retain=col,*$name \
    | sddsxref -pipe $twissFile -leave=* -transfer=param,* \
    | sddsprocess -pipe \
    "{-define=col,phi,x$name y$name atan2}" \
    "{-define=col,theta,x$name sqr betax / betax0 / y$name sqr betay / betay0 / + sqrt}" \
    "{-define=col,CS,1 theta 2 / tan sqr /}" \
    | sddssort -pipe -column=phi,incr -unique \
    | sddsprocess -pipe \
    -process=CS,integral,%sInteg,functionOf=phi \
    -define=col,pCentral0,pCentral \
    | sddscollapse -pipe \
    | sddsbreak -pipe -decreaseOf=s \
    | sddsprocess -pipe \
    -process=CSInteg,integ,%sInteg,functionOf=s \
    -proces=s,max,%sMax \
    -process=pCentral,first,%s0 \
    "{-define=param,c1,re_mks sqr c_mks * $P * sMax / pCentral0 sqr / kb_mks / 298 /}" \
    "{-define=param,decayRate,c1 CSIntegInteg * $fSZ2Sum *,units=1/s}" \
    "{-define=param,elasticScatteringLifetime,decayRate rec 3600 /,units=h}" \
    | sddscollapse -pipe \
    | sddsprocess -pipe=in $output \
    -process=elasticScatteringLifetime,ave,%sAve \
    -process=elasticScatteringLifetime,stand,%sStDev




