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

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)]

set usage "usage: longitCalcs -twiss <filename>"
set twiss ""
set args $argv
if {[APSStrictParseArguments {twiss}] || ![string length $twiss]} {
    puts stderr $usage
    exit 1
}
if ![file exists $twiss] {
    puts stderr "Not found: $twiss"
    exit 1
}

set dataItemList [list U0 Sdelta0 ex0 alphac pCentral revFreq sMax]
if [catch {exec sddsconvert $twiss -pipe=out -topage=1 \
             | sddsprocess -pipe -process=s,max,sMax \
             "-define=parameter,revFreq,sMax pCentral beta.p c_mks * / rec" \
             | sdds2stream -pipe \
             -parameter=[join $dataItemList ,]} dataList] {
    puts stderr "$dataList"
    exit 1
}
set index 0
foreach item $dataItemList {
    set ${item} [lindex $dataList $index]
    incr index
}
set revFreq1 [expr $revFreq]
set sMax1 [expr $sMax]

APSApplication . -name longitCalcs 
set harmonic 1296
set extradE 0
set harmonic2 3
set voltage2 0
set superPeriods 1
# initial start with overvoltage of 1.5
# to avoid problem with arcsin function
set voltage [expr $superPeriods*($U0+$extradE) * 1.5]
if {$voltage < 9} {set voltage 9}

APSLabeledEntry .voltage -parent .userFrame \
  -label "Voltage (MV): " -textVariable voltage -width 22
APSLabeledEntry .harmonic -parent .userFrame \
  -label "Harmonic number: " -textVariable harmonic -width 22 -type integer
APSLabeledEntry .sdelta -parent .userFrame \
    -label "Sigma P/P0: " -textVariable Sdelta0 -width 22
APSLabeledEntry .extradE -parent .userFrame \
    -label "Additional dE (MeV/period/turn): " -textVariable extradE -width 22
APSLabeledEntry .super -parent .userFrame \
    -label "Superperiods: " -textVariable superPeriods -width 22 
foreach item {voltage harmonic sdelta extradE super} {
    bind .userFrame.$item.entry <Return> compute
}

APSFrame .hharm -parent .userFrame  -label "Harmonic cavity" -width 22
set harmonicMode None
APSRadioButtonFrame  .mode -parent .userFrame.hharm.frame -label "Mode: " \
    -orientation horizontal -variable harmonicMode \
    -buttonList "None Lengthen Shorten" \
    -valueList "None Lengthen Shorten" \
    -commandList [list compute compute compute]
set harmonicRatio 4
set voltageRatio 0
set harmonicPhase 0
APSLabeledEntry .harmonicRatio -parent .userFrame.hharm.frame \
    -label "Freq. ratio" -textVariable harmonicRatio -type integer -width 22
APSLabeledEntry .voltageRatio -parent .userFrame.hharm.frame \
    -label "Volt. ratio" -textVariable voltageRatio -width 22 -type real
APSLabeledEntry .harmonicCavityPhase -parent .userFrame.hharm.frame \
    -label "Phase (deg)" -textVariable harmonicCavityPhase -width 22 -type real
bind .userFrame.hharm.frame.harmonicRatio.entry <Return> compute
bind .userFrame.hharm.frame.voltageRatio.entry <Return> compute
set lastHarmonicMode None

set synchTune 0
set bunchLength 0
set bunchDuration 0
set synchPhase 0
set synchPhaseDeg 0
set overVoltage 0
set rfAcceptance 0
set rfFrequency 0
set synchFreq 0
set coupling 0
APSLabeledOutput .length -parent .userFrame \
    -label "Circumference (m): " -textVariable sMax1 -width 22
APSLabeledOutput .revFreq -parent .userFrame \
    -label "Rev. frequency (MHz): " -textVariable revFreq1 -width 22
APSLabeledOutput .rfFreq -parent .userFrame \
    -label "RF frequency (MHz): " -textVariable rfFrequency -width 22
APSLabeledOutput .synchphase -parent .userFrame \
  -label "Synchrotron Phase: " -textVariable synchPhaseDeg -width 22
APSLabeledOutput .overvolt -parent .userFrame \
    -label "Over Voltage: " -textVariable overVoltage  -width 22
APSLabeledOutput .accept -parent .userFrame \
    -label "RF Acceptance: " -textVariable rfAcceptance -width 22
APSLabeledOutput .nus -parent .userFrame \
  -label "Synchrotron Tune: " -textVariable synchTune  -width 22
APSLabeledOutput .nuf -parent .userFrame \
  -label "Synchrotron Frequency (Hz): " -textVariable synchFreq  -width 22
APSLabeledOutput .blen -parent .userFrame \
  -label "Bunch length (m): " -textVariable bunchLength -width 22
APSLabeledOutput .bdur -parent .userFrame \
  -label "Bunch duration (s): " -textVariable bunchDuration -width 22
APSLabeledOutput .coupling -parent .userFrame \
  -label "Longitudinal coupling: " -textVariable coupling -width 22

proc compute {} {
    global dataItemList
    eval global $dataItemList bunchLength harmonic voltage synchPhaseDeg
    global overVoltage rfAcceptance rfFrequency synchFreq coupling extradE synchTune
    global bunchDuration superPeriods revFreq1 sMax1 sMax
    global harmonicRatio voltageRatio harmonicMode lastHarmonicMode harmonicCavityPhase
    global pi cMKS U1 sMax1 revFreq1 energy omegaRF rfFrequency

    if [expr $harmonicRatio<=1] {
        set harmonicMode None
        set harmonicRatio 1
    }

    set pi 3.141592653589793
    set cMKS 2.997924580000000e+08
    set U1 [expr $superPeriods*($U0+$extradE)]
    set sMax1 [expr $superPeriods*$sMax]
    set revFreq1 [expr $revFreq/$superPeriods]
    set energy [expr sqrt($pCentral*$pCentral+1)*0.51099906]
    set omegaRF [expr $harmonic*2*$pi*$revFreq1]
    set rfFrequency [expr $omegaRF/(2*$pi*1.0e6)]

    set overVoltage [expr $voltage/$U1]
    switch $harmonicMode {
        None {
            set voltageRatio 0
            set harmonicCavityPhase 0
            set synchPhase [expr asin($U1/$voltage)]
            set rfAcceptance [expr sqrt(2*$U1/($pi*$alphac*$harmonic*$energy)*(sqrt($overVoltage*$overVoltage-1)-acos(1/$overVoltage)))]
            set synchTune [expr sqrt($alphac*$harmonic*cos($synchPhase)*$voltage/(2*$pi*$energy))]
            set synchFreq [expr $synchTune*$revFreq1]
            set bunchLength [expr $Sdelta0*$cMKS*sqrt($alphac*$energy/($revFreq1*$omegaRF*cos($synchPhase)*$voltage))]
            set coupling [expr sqrt($alphac*$voltage*cos($synchPhase)*$omegaRF/($energy*$revFreq1*4.0))]
        }
        Lengthen {
            set m $harmonicRatio*1.0
            set f1 [expr $m*$m/($m*$m-1)/$overVoltage]
            set synchPhase [expr asin($f1)]
            set k [expr sqrt(1/($m*$m) - 1/(($m*$m-1)*pow($overVoltage,2)))]
            set voltageRatio $k
            set harmonicPhase [expr asin(-1/($overVoltage*($m*$m-1)*$k))]
            set harmonicCavityPhase [expr $harmonicPhase*180/$pi]
            if [expr $harmonicCavityPhase<0] {
                set harmonicCavityPhase [expr $harmonicCavityPhase+360]
            }
            #set check [expr $voltage*(sin($synchPhase) + $voltageRatio*sin($harmonicPhase))-$U1]
            #puts stderr "check: $check"
            set synchTune0 [expr sqrt($alphac*$harmonic*cos($synchPhase)*$voltage/(2*$pi*$energy))]
            set synchFreq ?
            set synchTune ?
            set rfAcceptance ?
            set coupling ?
            set bunchLength [expr 2/3.6256*pow(3./($m*$m-1), 0.25)*sqrt(2*$harmonic*$alphac*$Sdelta0/$synchTune0)*$cMKS/$omegaRF]
        }
        Shorten {
            if [expr $voltageRatio<0] {
                set voltageRatio 0
            }
            set synchPhase [expr asin($U1/$voltage)]
            set synchTune0 [expr sqrt($alphac*$harmonic*cos($synchPhase)*$voltage/(2*$pi*$energy))]
            set synchTune [expr $synchTune0*sqrt(1 + $harmonicRatio*$voltageRatio/cos($synchPhase))]
            set synchFreq [expr $synchTune*$revFreq1]
            set bunchLength [expr $Sdelta0*$cMKS*sqrt($alphac*$energy/($revFreq1*$omegaRF*cos($synchPhase)*$voltage))/sqrt(1 + $harmonicRatio*$voltageRatio/cos($synchPhase))]
            set coupling ?
        }
    }

    set synchPhaseDeg [expr 180-180.0*$synchPhase/$pi]
    set lastHarmonicMode $harmonicMode
    set bunchDuration [expr $bunchLength/$cMKS]

    computeBunchShape

    update
}

proc computeBunchShape {} {
    global dataItemList
    eval global $dataItemList bunchLength harmonic voltage synchPhaseDeg twiss
    global overVoltage rfAcceptance rfFrequency synchFreq coupling extradE synchTune energy
    global bunchDuration superPeriods revFreq1 sMax1 sMax
    global harmonicRatio voltageRatio harmonicMode lastHarmonicMode harmonicCavityPhase
    global pi cMKS U1 sMax1 revFreq1 energy omegaRF rfFrequency

    exec sddssequence -pipe=out -define=t,type=double,units=s -sequence=begin=[expr -0.25e-6/$rfFrequency],end=[expr 0.25e-6/$rfFrequency],n=10000 \
        | sddsprocess -pipe=in $twiss.longit \
        "-define=param,C1,[expr $alphac*$cMKS*$voltage/($energy*$sMax1*$omegaRF)]" \
        "-define=param,alphac,$alphac" \
        "-define=param,Sdelta0,$Sdelta0" \
        "-define=param,phis,$synchPhaseDeg dtor" \
        "-define=param,k,$voltageRatio" \
        "-define=param,n,$harmonicRatio,type=short" \
        "-define=param,phih,$harmonicCavityPhase dtor" \
        "-define=param,omega,$omegaRF" \
        "-define=col,phi,omega t *" \
        "-define=col,Phi1,phis cos  phi     phis + cos -" \
        "-define=col,Phi2,phih cos  phi n * phih + cos - k * n /" \
        "-define=col,Phi3,phis sin  phih sin k * + phi *" \
        "-define=col,Phi,Phi1 Phi2 + Phi3 - C1 *" \
        "-define=col,rhoArg,Phi alphac sqr / Sdelta0 sqr / chs" \
        -process=rhoArg,min,%sMin \
        "-define=col,rho,rhoArg chs exp" \
        -process=rho,integ,rhoInteg,functionOf=t \
        "-redefine=col,rho,rho rhoInteg /" \
        -process=t,rms,tRms,weight=rho

    exec sddsplot -column=t,rho $twiss.longit -title=@tRms &
}

compute
APSButton .compute -parent .userFrame -text Compute -command compute
