#!/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 auto_path [linsert $auto_path 0 /usr/local/oag/apps/lib/$env(HOST_ARCH)]

set usage "usage: longitCalcs \[-gui 0\] \[-output <filename>\] -twiss <filename> \[-voltage <MV> -harmonic <number> -superPeriods <number>(1)\]"
set twiss  /home/helios/oagData/sr/lattices/default/aps.twi
set voltage 4.0
set harmonic 1296
set gui 1
set output ""
set superPeriods 1
set dEExtra 0.0
set beamCurrent 0.0
set args $argv
if {[APSStrictParseArguments {output superPeriods twiss voltage harmonic gui dEExtra beamCurrent}] || ![string length $twiss]} {
    puts stderr $usage
    exit 1
}

set beamCurrentPV S-DCCT:CurrentM
set rfVoltagePV SR:C:cavTotalGapVolt
set rfVoltagePVFactor 1e-3
set hrfVoltagePV S38-BLS:LLRF:CavitySampleCalMagM
set hrfVoltagePVFactor 1.0
set IDLossPV ID:TotalEnergyLossM
set IDLossPVFactor 1.0
set hVoltage 0
set hh 4
set hcOpt 0
set hphase 360

if $gui==0 {
    wm withdraw .
    if ![string length $output] {
        puts stderr "Supply output file with -output option"
        exit 1
    }
}
if [file exists $output] {
    puts stderr "In use: $output"
    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]
set gamma [expr sqrt(1+$pCentral*$pCentral)]
set eta [expr $alphac-1./pow($gamma,2)]

set extradE $dEExtra
set harmonic2 $hh
set voltage2 $hVoltage
if [expr $voltage<0] {
    # initial start with overvoltage of 1.5
    # to avoid problem with arcsin function
    set voltage [expr $superPeriods*($U0+$extradE) * 1.5]
}

set harmonicRatio $harmonic2 
set harmonicVoltage $hVoltage
set voltageRatio 0
set harmonicMode None
if [expr $harmonicVoltage!=0] {
    set harmonicMode Custom
}
if $hcOpt {
    set harmonicMode OptLengthen
}
if  [string length $hrfVoltagePV] {
    set harmonicMode PV
}
set harmonicPhase $hphase
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
set fractionalSpan 0.25
set harmonicPhaseDeg $hphase
set harmonicROverQ 103.0
set harmonicLoadedQ 7.3e6
set harmonicDeltaFreq 15.0
set usePVValues 0
set nCavities 12
set QLoaded 12000
set ROverQa 228.0

set statusText ""
proc setStatus {text} {
    global statusText
    set statusText "$text"
    update
}

if $gui {
    APSApplication . -name longitCalcs 
    APSScrolledStatus .status -parent .userFrame -width 50  -textVariable statusText
    APSLabeledEntry .current -parent .userFrame \
        -label "Beam current (mA): " -textVariable beamCurrent -width 22
    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 
    APSLabeledEntry .nCavities -parent .userFrame \
        -label "Number of cavities: " -textVariable nCavities -width 22 
    APSLabeledEntry .loadedQ -parent .userFrame \
        -label "Loaded Q: " -textVariable QLoaded -width 22 
    APSLabeledEntry .rOverQMain -parent .userFrame \
        -label "R/Q (Ohm): " -textVariable ROverQa -width 22 
    foreach item {voltage harmonic sdelta extradE super} {
        bind .userFrame.$item.entry <Return> compute
    }

    if {[string length $IDLossPV] || [string length $rfVoltagePV]} {
        set usePVValues 1
        APSRadioButtonFrame .rbpv -parent .userFrame -label "Use PV values" \
          -orientation horizontal -variable usePVValues -buttonList "Yes No" \
          -valueList "1 0" -commandList [list compute compute] \
          -contextHelp "If checked, use PV values for main rf voltage and extra energy loss"
    }           

    APSFrame .hharm -parent .userFrame  -label "Harmonic cavity" -width 22
    if [string length $hrfVoltagePV] {
        APSRadioButtonFrame  .mode -parent .userFrame.hharm.frame -label "Mode: " \
          -orientation horizontal -variable harmonicMode \
          -buttonList [list None "Opt.Lengthen" "Custom" "PV"] \
          -valueList [list None OptLengthen Custom PV] \
          -commandList [list compute compute compute compute]
    } else {
        APSRadioButtonFrame  .mode -parent .userFrame.hharm.frame -label "Mode: " \
          -orientation horizontal -variable harmonicMode \
          -buttonList [list None "Opt.Lengthen" "Custom"] \
          -valueList [list None OptLengthen Custom] \
          -commandList [list compute compute compute]
    }
    
    APSLabeledEntry .hvoltage -parent .userFrame.hharm.frame \
        -label "Voltage (kV)" -textVariable harmonicVoltage -width 22 -type real
    APSLabeledEntry .harmonicPhaseDeg -parent .userFrame.hharm.frame \
        -label "Phase (deg)" -textVariable harmonicPhaseDeg -width 22 -type real
    APSLabeledEntry .roverq -parent .userFrame.hharm.frame \
        -label "R/Q (Ohm)" -textVariable harmonicROverQ -width 22 -type real
    APSLabeledEntry .ql -parent .userFrame.hharm.frame \
        -label "Loaded Q" -textVariable harmonicLoadedQ -width 22 -type real
    APSLabeledEntry .df -parent .userFrame.hharm.frame \
        -label "Freq. Offset (kHz)" -textVariable harmonicDeltaFreq -width 22 -type real
    bind .userFrame.hharm.frame.hvoltage.entry <Return> "compute -harmonicVoltageUpdated 1"
    bind .userFrame.hharm.frame.harmonicPhaseDeg.entry <Return> compute
    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 "Synchronous Phase (deg): " -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 .robinsonThreshold -parent .userFrame \
        -label "Robinson threshold (mA): " -textVariable RobinsonThreshold -width 22
    #APSLabeledOutput .coupling -parent .userFrame \
    #    -label "Longitudinal coupling: " -textVariable coupling -width 22
    APSLabeledEntry .span -parent .userFrame \
        -label "Calculation half-span: " -textVariable fractionalSpan -width 22 \
        -contextHelp "Half-span of the bunch distribution calculation in units of the rf period. Changing this may resolve numerical problems for long bunches."
} 

set outputIndex -1
proc compute {args} {
    set plot 0
    set noPVUpdates 0
    set harmonicVoltageUpdated 0
    APSStrictParseArguments {plot harmonicVoltageUpdated noPVUpdates}
    global dataItemList gui output outputIndex
    eval global $dataItemList bunchLength harmonic voltage synchPhaseDeg eta gamma synchPhase
    global overVoltage rfAcceptance rfFrequency synchFreq coupling extradE synchTune
    global bunchDuration superPeriods revFreq1 sMax1 sMax
    global harmonicRatio harmonicMode lastHarmonicMode harmonicPhaseDeg harmonicVoltage voltageRatio
    global harmonicROverQ harmonicLoadedQ harmonicDeltaFreq
    global pi cMKS U1 sMax1 revFreq1 energy omegaRF rfFrequency output
    global rfVoltagePV rfVoltagePVFactor IDLossPV IDLossPVFactor usePVValues hrfVoltagePV hrfVoltagePVFactor
    global beamCurrentPV beamCurrent RobinsonThreshold

    if {$usePVValues && !$noPVUpdates} {
        if [string length $beamCurrentPV] {
            set beamCurrent [exec cavget -list=$beamCurrentPV]
        }
        if [string length $rfVoltagePV] {
            set voltage [expr [exec cavget -list=$rfVoltagePV]*$rfVoltagePVFactor]
        }
        if [string length $IDLossPV] {
            set extradE [expr [exec cavget -list=$IDLossPV]*$IDLossPVFactor]
        }
        if {[string length $hrfVoltagePV] && [string compare $harmonicMode "PV"]==0} {
            set harmonicVoltage [expr [exec cavget -list=$hrfVoltagePV]*$rfVoltagePVFactor*1e6]
            set harmonicVoltageUpdated 1
        }
        setStatus "Updating values from PVs"
    }

    if [expr $harmonicRatio<=1] {
        set harmonicMode None
        set harmonicRatio 1
    }
    if [expr $voltage<=0] {
        setStatus "Rf voltage is invalid"
        return
    }
    set voltageRatio [expr $harmonicVoltage/$voltage/1e3]
    
    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 $gamma*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 harmonicPhase 0.0
            set harmonicPhaseDeg 0
            set synchPhase [expr $pi-asin($U1/$voltage)]
            set rfAcceptance [expr sqrt(2*$U1/($pi*$eta*$harmonic*$energy)*(sqrt($overVoltage*$overVoltage-1)-acos(1/$overVoltage)))]
            set synchTune [expr sqrt(-$eta*$harmonic*cos($synchPhase)*$voltage/(2*$pi*$energy))]
            set synchFreq [expr $synchTune*$revFreq1]
            set bunchLength [expr $Sdelta0*$cMKS*sqrt($eta*$energy/(-$revFreq1*$omegaRF*cos($synchPhase)*$voltage))]
            set coupling [expr sqrt(-$eta*$voltage*cos($synchPhase)*$omegaRF/($energy*$revFreq1*4.0))]
            set bunchDuration [expr $bunchLength/$cMKS]
        }
        OptLengthen {
            set m $harmonicRatio*1.0
            set f1 [expr $m*$m/($m*$m-1)/$overVoltage]
            set synchPhase [expr $pi-asin($f1)]
            #puts stderr "cos(synchPhase) = [expr cos($synchPhase)]"
            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 harmonicPhaseDeg [expr $harmonicPhase*180/$pi]
            if [expr $harmonicPhaseDeg<0] {
                set harmonicPhaseDeg [expr $harmonicPhaseDeg+360]
            }
            #set check [expr $voltage*(sin($synchPhase) + $voltageRatio*sin($harmonicPhase))-$U1]
            #puts stderr "check: $check"
            set synchTune0 [expr sqrt(-$eta*$harmonic*cos($synchPhase)*$voltage/(2*$pi*$energy))]
            set bunchLength [expr 2/3.6256*pow(3./($m*$m-1), 0.25)*sqrt(2*$harmonic*$eta*$Sdelta0/$synchTune0)*$cMKS/$omegaRF]
            set bunchDuration [expr $bunchLength/$cMKS]
            set synchFreq [expr 0.128*$eta*$Sdelta0/$bunchDuration]
            set synchTune [expr $synchFreq/$revFreq1]
            set rfAcceptance ?
            set coupling ?
        }
        PV -
        Custom {
            if [expr $voltageRatio<0] {
                set voltageRatio 0
            }
            set harmonicPhase [expr $harmonicPhaseDeg*$pi/180.]
            set U2 [expr $U1-$voltageRatio*$voltage*sin($harmonicPhase)]
            #setStatus "U2 = $U2"
            set synchPhase [expr $pi-asin($U2/$voltage)]
            #set check [expr $voltage*(sin($synchPhase) + $voltageRatio*sin($harmonicPhase))-$U1]
            #puts stderr "check: $check"
            set synchPhaseDeg [expr $synchPhase*180/$pi]
            set synchTune0 [expr sqrt(-$eta*$harmonic*cos($synchPhase)*$voltage/(2*$pi*$energy))]
            set synchTune ?
            set synchFreq ?
            set bunchLength [expr $Sdelta0*$cMKS*sqrt($eta*$energy/(-$revFreq1*$omegaRF*cos($synchPhase)*$voltage))/sqrt(1 + $harmonicRatio*$voltageRatio/(-cos($synchPhase)))]
            set coupling ?
        }
    }
    set harmonicVoltage [expr $voltageRatio*$voltage*1e3]

    set synchPhaseDeg [expr 180.0*$synchPhase/$pi]
    set lastHarmonicMode $harmonicMode

    computeBunchShape -plot $plot
    findFixedPoints
    computeBucketHH
    computeRobinsonThreshold

    if $plot {
        set Trf [expr 1e-6/$rfFrequency]
        exec sddssequence -pipe=out -define=t,units=s -sequence=begin=[expr -$Trf*0.75],end=[expr $Trf*0.75],n=1000 \
          | sddsprocess -pipe=in $output.volt  \
            "-define=col,mainPhase,$omegaRF t * $synchPhase +" \
            "-define=col,harmPhase,$omegaRF t * $harmonicRatio * $harmonicPhase +" \
            "-define=col,VMain,mainPhase sin $voltage *,units=MV" \
            "-define=col,VHarm,harmPhase sin $voltage * $voltageRatio *,units=MV" \
            "-define=col,VTotal,VMain VHarm +,units=MV"  
        set ts [expr $synchPhase/(2*$pi)*$Trf*1e9]
        exec sddsplot -convert=col,t,ns,s,1e9 -graph=line,vary,type=1,thick=2 "-ylabel=Voltage (MV)" \
          -column=t,VMain -legend=spec=Main $output.volt \
          -column=t,VHarm -legend=spec=Harm. $output.volt \
          -column=t,VTotal -legend=spec=Total $output.volt \
          -axes -drawline=p0v=0.4,p1v=0.6,y0v=[expr $U0+$extradE],y1v=[expr $U0+$extradE],linetype=5 &
        # -drawline=linetype=4,q0v=0,q1v=1,x0v=$ts,x1v=$ts 
    }
    update
    if [string length $output] { 
        if [string compare $synchTune ?]==0 {
            set synchTuneNum -1
            set synchFreqNum -1
        } else {
            set synchTuneNum $synchTune
            set synchFreqNum $synchFreq
        }
        exec sddsmakedataset $output.tmp \
            -column=MainVoltage,units=MV,type=double -data=$voltage \
            -column=MainHarmonic,type=long -data=$harmonic \
            -column=Sdelta0,type=double -data=$Sdelta0 \
            -column=U0,type=double,units=MeV -data=$U0 \
            -column=dU,type=double,units=MeV -data=$extradE \
            -column=Superperiods,type=long -data=$superPeriods \
            -column=HHRatio,type=short -data=$harmonicRatio \
            -column=VhRatio,type=double -data=$voltageRatio \
            -column=Vh,type=double,units=MV -data=[expr $voltageRatio*$voltage] \
            -column=HPhase,type=double,units=deg -data=$harmonicPhaseDeg \
            -column=Circumference,type=double,units=m -data=$sMax1 \
            -column=RfFrequency,type=double,units=MHz -data=$rfFrequency \
            -column=MainPhase,type=double,units=deg -data=$synchPhaseDeg \
            -column=SynchrotronFrequency,type=double,units=Hz -data=$synchFreqNum \
            -column=SynchrotronTune,type=double -data=$synchTuneNum \
            -column=RfAcceptance,type=double -data=[expr $rfAcceptance/100.0] \
            -column=BunchLength,type=double,units=m -data=$bunchLength \
            -column=BunchDuration,type=double,units=s -data=$bunchDuration \
            -column=BeamCurrent,type=double,units=mA -data=$beamCurrent \
            -column=RobinsonThreshold,type=double,units=mA -data=$RobinsonThreshold 
        if [file exists $output] {
            exec sddscombine $output $output.tmp $output.tmp2 -merge -over
            file rename -force $output.tmp2 $output
            file delete $output.tmp
        } else {
            file rename $output.tmp $output
        }
    }
}

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

    #puts stderr "findFixedPoints: $voltage, $voltageRatio $harmonicRatio $synchPhaseDeg $harmonicPhaseDeg $U1"
    set dataList [exec sddssequence -pipe=out -define=dphi -sequence=begin=[expr -2*$pi],end=[expr 2*$pi],n=1000000 \
		      | sddsprocess -pipe -clip=1,0 \
		      "-define=param,V1,$voltage" "-define=param,Vh,$voltage $voltageRatio *" \
		      "-define=param,h,$harmonicRatio" \
		      "-define=param,phi1,$synchPhaseDeg pi * 180 /" "-define=param,phih,$harmonicPhaseDeg pi * 180 /" \
		      "-define=param,U0,$U1" \
		      "-define=column,error,U0 V1 phi1 dphi + sin * - Vh phih dphi h * + sin * -" \
		      "-define=column,deriv,V1 phi1 dphi + cos * Vh phih dphi h * + cos * h * + chs" \
		      | sddszerofind -pipe -zeroesOf=error -columns=dphi,deriv \
		      | sddsprocess -pipe "-test=column,deriv 0 <" \
		      | sdds2stream -pipe -column=dphi]
    global fixedPoint
    if [llength $dataList]>1 {
	set fixedPoint(1) [lindex $dataList 0]
	set fixedPoint(2) [lindex $dataList 1]
	set fixedPoint(valid) 1
    } else {
	set fixedPoint(valid) 0
    }
}

proc computeBucketHH {args} {
    global dataItemList gui
    eval global $dataItemList bunchLength harmonic voltage synchPhaseDeg twiss eta gamma
    global overVoltage rfAcceptance rfFrequency synchFreq coupling extradE synchTune energy
    global bunchDuration superPeriods revFreq1 sMax1 sMax
    global harmonicRatio voltageRatio harmonicMode lastHarmonicMode harmonicPhaseDeg
    global pi cMKS U1 sMax1 revFreq1 energy omegaRF rfFrequency
    global fixedPoint rfAcceptance

    if $fixedPoint(valid)==0 return
    set fp $fixedPoint(1)
    set phi1 [expr $synchPhaseDeg*$pi/180]
    set phih [expr $harmonicPhaseDeg*$pi/180]
    set Vh [expr $voltageRatio*$voltage]
    set W1 [expr $voltage*(cos($phi1) - cos($phi1+$fp))]
    set W2 [expr $Vh*(cos($phih) - cos($phih+$fp*$harmonicRatio))/$harmonicRatio]
    set W3 [expr $U1*$fp]
    set W [expr $W1+$W2-$W3]
    if [expr $W>=0] {
        set rfAcceptance 0
    } else {
	set rfAcceptance [expr sqrt(-$W/($pi*$harmonic*$eta*$energy))*100]
    }
    
}

proc computeBunchShape {args} {
    set plot 0
    APSStrictParseArguments {plot}
    global dataItemList
    eval global $dataItemList bunchLength harmonic voltage synchPhaseDeg twiss eta gamma env
    global overVoltage rfAcceptance rfFrequency synchFreq coupling extradE synchTune energy
    global bunchDuration superPeriods revFreq1 sMax1 sMax
    global harmonicRatio voltageRatio harmonicMode lastHarmonicMode harmonicPhaseDeg
    global harmonicROverQ harmonicLoadedQ harmonicDeltaFreq
    global pi cMKS U1 sMax1 revFreq1 energy omegaRF rfFrequency fractionalSpan

    set outFile  $twiss.longit
    if [catch {exec touch $outFile} result] {
        set outFile [file tail $twiss].longit
        if [catch {exec touch $outFile} result] {
            set outFile /tmp/$env(USER)-[file tail $twiss].longit
        }
    }
    exec sddssequence -pipe=out -define=t,type=double,units=s -sequence=begin=[expr -$fractionalSpan*1e-6/$rfFrequency],end=[expr $fractionalSpan*1e-6/$rfFrequency],n=10000 \
        | sddsprocess -pipe=in $outFile \
        "-define=param,C1,[expr $eta*$cMKS*$voltage/($energy*$sMax1*$omegaRF)]" \
        "-define=param,alphac,$alphac" \
        "-define=param,eta,$eta" \
        "-define=param,Sdelta0,$Sdelta0" \
        "-define=param,phis,$synchPhaseDeg dtor" \
        "-define=param,k,$voltageRatio" \
        "-define=param,n,$harmonicRatio,type=short" \
        "-define=param,phih,$harmonicPhaseDeg 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 eta 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,ave,tAve,weight=rho \
        -process=t,rms,tRms,weight=rho \
        "-define=col,rhoSqr,rho sqr" \
        -process=rhoSqr,integ,rhoSqrInteg,functionOf=t \
        "-define=col,FKReal,rho omega n * t * cos *" \
        "-define=col,FKImag,rho omega n * t * sin * chs" \
        "-process=FKReal,integ,FReal,functionOf=t" \
        "-process=FKImag,integ,FImag,functionOf=t" \
        "-define=parameter,FAbs,FReal sqr FImag sqr + sqrt"     

    set bunchDuration [exec sdds2stream -parameter=tRms $outFile]
    set bunchOffset [exec sdds2stream -parameter=tAve $outFile]
    set bunchLength [expr $bunchDuration*$cMKS]
    set F [exec sdds2stream -parameter=FAbs $outFile]

    if $plot {
        exec sddsplot -column=t,rho $outFile "-title=@tRms,edit=i/Rms duration (s): /" \
            "-topline=@rhoSqrInteg,edit=i/Integral sqr(rho): /" &
    }
    return [list $bunchDuration $bunchOffset $F]
}

proc outputForElegantDialog {} {
    global gui
    global code outputFile rfMainName rfHarmonicName harmonicMode eta gamma
    global voltage rfFrequency synchPhaseDeg voltageRatio harmonicRatio harmonicPhaseDeg
    if $gui {
        APSDialogBox .output -modal 0 -okCommand "set code 1" -cancelCommand "set code 0" \
            -name "Elegant output"
        APSLabeledEntry .le1 -parent .output.userFrame -label "Filename: " \
            -textVariable outputFile -width 80 
        APSLabeledEntry .le2 -parent .output.userFrame -label "Main rf element name: " \
            -textVariable rfMainName -width 80
        if {[string compare $harmonicMode None] && $harmonicRatio>1} {
            APSLabeledEntry .le3 -parent .output.userFrame -label "Harmonic rf element name: " \
                -textVariable rfHarmonicName -width 80 
            set doHarmonic 1
        } else {
            set doHarmonic 0
            set rfHarmonicName ""
        }
        update
        tkwait variable code
        if {$code && [string length $outputFile] && [string length $rfMainName] && (!$doHarmonic || [string length $rfHarmonicName])} {
            compute -plot 0
            lappend ElementNameList $rfMainName $rfMainName $rfMainName 
            lappend ElementParameterList VOLT FREQ PHASE
            lappend ParameterValueList [expr $voltage*1e6] [expr $rfFrequency*1e6] $synchPhaseDeg
            if $doHarmonic {
                lappend ElementNameList $rfHarmonicName $rfHarmonicName $rfHarmonicName
                lappend ElementParameterList VOLT FREQ PHASE
                lappend ParameterValueList [expr $voltage*1e6*$voltageRatio] [expr $rfFrequency*1e6*$harmonicRatio] $harmonicPhaseDeg
            }
            exec sddsmakedataset $outputFile \
                -column=ElementName,type=string -data=[join $ElementNameList ,] \
                -column=ElementParameter,type=string -data=[join $ElementParameterList ,] \
                -column=ParameterValue,type=double -data=[join $ParameterValueList ,]
        }
    }
}

proc iterateBunchLength {} {
    global dataItemList gui output outputIndex
    eval global $dataItemList bunchLength harmonic voltage synchPhaseDeg eta gamma
    global overVoltage rfAcceptance rfFrequency synchFreq coupling extradE synchTune
    global bunchDuration superPeriods revFreq1 sMax1 sMax
    global harmonicRatio harmonicMode lastHarmonicMode harmonicPhaseDeg harmonicVoltage voltageRatio
    global harmonicROverQ harmonicLoadedQ harmonicDeltaFreq
    global pi cMKS U1 sMax1 revFreq1 energy omegaRF rfFrequency output
    global rfVoltagePV rfVoltagePVFactor IDLossPV IDLossPVFactor usePVValues hrfVoltagePV hrfVoltagePVFactor
    global beamCurrentPV beamCurrent

    setStatus "Performing iteration"

    compute -plot 0 

    set St $bunchDuration
    set Ct 0.0
    set lastSt 0
    set lastCt 0

    set left 100
    set cf 0.25
    set Vh [expr $harmonicVoltage*1e3]
    set I [expr $beamCurrent/1e3]
    set RoQh $harmonicROverQ
    set Qlh $harmonicLoadedQ
    set nh 1
    set hh $harmonicRatio
    set Vm [expr $voltage*1e6]
    set frf [expr $rfFrequency*1e6]

    set cycles 0
    set F 1.0
    while {([expr (abs($St-$lastSt)/$St)>1e-4] && $left>0) || $cycles<2} {
        # Compute detuning of HHC
        #set psih [expr acos($Vh/($nh*$I*$RoQh*$Qlh*exp(-pow($omegaRF*$hh*$St, 2)/2)))]
        set psih [expr acos($Vh/($nh*$I*$RoQh*$Qlh*$F))]
        set dfh  [expr tan($psih)*$frf*$hh/(2*$Qlh)]
        set harmonicCavityPhase [expr $psih-$pi/2.+$Ct*$frf*$hh*2*$pi]
        set harmonicPhaseDeg [expr 360+$harmonicCavityPhase*180/$pi]
        set harmonicDeltaFreq [expr $dfh/1e3]

        # Compute synchronous phase for the main cavity
        set phis [expr $pi-asin((1e6*$U1-$Vh*sin($harmonicCavityPhase))/$Vm)]
        set synchPhaseDeg [expr $phis*180/$pi]
        #setStatus "synchPhaseDeg = $synchPhaseDeg"

        # Compute the bunch shape and offset
        set lastSt $St
        set lastCt $Ct
        set dataList [computeBunchShape]
        set St [lindex $dataList 0]
        set Ct [lindex $dataList 1]
        set St [expr ($cf*$St+(1-$cf)*$lastSt)]
        set Ct [expr ($cf*$Ct+(1-$cf)*$lastCt)]
        set F [lindex $dataList 2]
        incr left -1
        incr cycles
        setStatus "St=$St after $cycles iterations"
        if $cycles==50 {
            set $cf [expr $cf/2.0]
        }
    }
    compute -plot 0 -noPVUpdates 1
}

proc computeRobinsonThreshold {} {
    global dataItemList gui output outputIndex
    eval global $dataItemList bunchLength harmonic voltage synchPhaseDeg eta gamma synchPhase
    global overVoltage rfAcceptance rfFrequency synchFreq coupling extradE synchTune
    global bunchDuration superPeriods revFreq1 sMax1 sMax
    global harmonicRatio harmonicMode lastHarmonicMode harmonicPhaseDeg harmonicVoltage voltageRatio
    global harmonicROverQ harmonicLoadedQ harmonicDeltaFreq
    global pi cMKS U1 sMax1 revFreq1 energy omegaRF rfFrequency output
    global rfVoltagePV rfVoltagePVFactor IDLossPV IDLossPVFactor usePVValues hrfVoltagePV hrfVoltagePVFactor
    global beamCurrentPV beamCurrent RobinsonThreshold
    global nCavities QLoaded ROverQa

    set VCavity [expr (1e6*$voltage)/$nCavities]
    set Idc [expr $beamCurrent/1e3]
    set pi [expr 4*atan(1)]
    set ff [expr exp(-pow($rfFrequency*2e6*$pi*$bunchDuration,2))]
    #setStatus "ff: $ff  synchPhase: $synchPhase"
    for {set i 0} {$i<10} {incr i} {
        set phi_z_opt [expr atan($ff*$Idc*$ROverQa*$QLoaded*cos($synchPhase)/$VCavity)]
        set Idc [expr 2*$VCavity*cos($synchPhase)/sin(2*$phi_z_opt)/($ROverQa*$QLoaded)/$ff]
    }
    set RobinsonThreshold [expr $Idc*1e3]
}

compute
if $gui {
    APSButton .compute -parent .userFrame -text Compute -command "compute -plot 0"
    APSButton .computeSC -parent .userFrame -text "Compute Self Consistent" -command "iterateBunchLength"
    APSButton .plot -parent .userFrame -text Plot -command "compute -plot 1"
    APSButton .output  -parent .userFrame -text ">Elegant" -command outputForElegantDialog
} else {
    exit 0
}
