#
# 
#

DEBUG=0
DEBUG2=0

if test $# -eq 0 ; then 
   outfile=`echo $0 | sed 's/.mkin$/.in/g'`
else
   outfile=`echo $0 | sed 's/.mkin$/.in.deb/g'`
fi

if test $# -ge 1 ; then
   DEBUG=1
fi

if test $# -ge 2 ; then
   DEBUG2=1
fi

echo 1 |awk 'BEGIN{

PI = atan2(1.,1.)*4.0 ; 
ECHARGE = 1.6021892e-19;
c0=2.99792458e8 ;

DEBUG='$DEBUG' 
DEBUG2='$DEBUG2' 

if(DEBUG){print "#DEBUG,DEBUG2 = ", DEBUG,DEBUG2}

Harmonics = 2*2*3*7*29        #Harmonics = 2436

NPack = 2*3*7
NBunch = Harmonics/NPack
NBunchWithQ = NBunch

NParticle=60 ; 

NTurn = 80000 ; 
SamplingStep = 100 

Eref_eV = 8e9 ; 
FreqAcc = 508.579e6; 
frev = FreqAcc/Harmonics 
Tref = 1/frev

Freqref = 1.0/Tref
NuBeta = 51.22 ; 
Eta =  1.46e-4 ; 
Je = 2.0 

Chrom0 = 2 
Chrom1 = 0 
ChromFreq = 0
ChromPhase = 0 ;

U0  = 9.2e6 ; 
Uid = 0 ;

Emittance_Natural = 6.99e-9 
nde_rms_Natural = 1.094e-3
de_rms_Natural = nde_rms_Natural * Eref_eV ;

LDamping = 1 ;
LExcitation = 1 ;

TDamping = 0 ;
TExcitation = 0 ;

Beta=10; 
Alpha=0.0 ;

Current=0.1 ;

NLDevice = 0 ;

NTDevice = 0 ;

if(DEBUG) print "# NBunch        = ", NBunch
if(DEBUG) print "# NBunchWithQ   = ", NBunchWithQ
if(DEBUG) print "# NParticle     = ", NParticle
if(DEBUG) print "# NTurn         = ", NTurn

if(DEBUG) print "# Eref_eV, Tref = ", Eref_eV, Tref
if(DEBUG) print "# Current       = ", Current 
if(DEBUG) print "# NuBeta        = ", NuBeta
if(DEBUG) print "# Eta           = ", Eta
if(DEBUG) print "# Je            = ", Je 
if(DEBUG) print "# U, U0, Uid    = ", U0+Uid, U0, Uid
if(DEBUG) print "# Beta, Alpha   = ", Beta, Alpha
if(DEBUG) print "# Chrom0 Chrom1 ChromFreq = ",Chrom0,Chrom1,ChromFreq

#***************************************************8
#
# Acceleration Resonator 
#
#***************************************************8

     NAcc = 1  ;
 NAccReal = 24 ;
#    QAcc0 = 43909 ; 
#    RAcc0 = 3394589 * NAccReal ;
    QAcc0 = 40000 ;
    RAcc0 = 3e6 * NAccReal ;
 RperQAcc =  RAcc0/QAcc0 ; 
      Vc0 = 12e6 ;  
 PhaseAcc = 0
 BetaC    = 2.0

if(DEBUG) print "# NAcc, NAccReal = ", NAcc, NAccReal 
if(DEBUG) print "# QAcc0, RAcc0   = ", QAcc0, RAcc0 
if(DEBUG) print "# Vc0            = ", Vc0
if(DEBUG) print "# FreqAcc        = ", FreqAcc 
if(DEBUG) print "# BetaC          = ", BetaC

#***************************************************8
#
# ======== Get PhiS, TauS, NuS
#
#    Synchrotron Osicllation Parameters
#
#***************************************************8

U = U0 + Uid
OmegaAcc = 2*PI*FreqAcc
sinPhiS = U/Vc0
cosPhiS = sqrt( 1 - sinPhiS*sinPhiS)
PhiS = atan2(sinPhiS,cosPhiS)
OmegaS = sqrt( Eta*Vc0/(Eref_eV*Tref)*OmegaAcc*cosPhiS )
FreqS = OmegaS/(2*PI)
NuS = FreqS/Freqref
tauS = PhiS/OmegaAcc 

if(DEBUG) print "# NuS, PhiS(deg.) = ", NuS, PhiS/(2*PI)*360 

PhaseAcc2 = PhiS ;

#***************************************************8
#
# Longitudinal Resonator
#
#***************************************************8

NLRes = 1 ;

NLResReal = 32 ; 
# TM011 of Bell Shape
   RperQ_L = 32.35 * NLResReal ;
       Q_L = 2000  ;
   FreqRes_L_0 = 910e6 

      h_HOM_L = int( FreqRes_L_0 / frev ) 
        P_L = int( h_HOM_L / Harmonics )
       Mu_L = h_HOM_L - P_L*Harmonics
 FreqRes_L = (P_L*Harmonics + Mu_L + NuS) * frev

       R_L = RperQ_L * Q_L

 if( DEBUG ) print "# Longitudinal" 
 if(DEBUG) print "# FreqRes_L_0, FreqRes_L = ", FreqRes_L_0, FreqRes_L 
 if(DEBUG) print "# h_HOM_L, Mu_L, P_L = ", h_HOM_L, Mu_L, P_L ; 
 if(DEBUG) print "# R_L, Q_L = ", R_L, Q_L ; 

#***************************************************8
#
# Transverse Resonator
#
#***************************************************8

NTRes = 1 ;

NTResReal=1 ;

# TM111 of Bell Shape
        Q_T = 40000 ;
    RperQ_T = 313.25 * NTResReal ;
   FreqRes_T_0 = 1083.0e6

       h_HOM_T = int( FreqRes_T_0 / frev ) 
        P_T = int( h_HOM_T / Harmonics ) + 1 
       Mu_T =  P_T * Harmonics - h_HOM_T - int( NuBeta )
  FreqRes_T = ( P_T * Harmonics - Mu_T - NuBeta ) * frev

        R_T = RperQ_T * Q_T

 if( DEBUG ) print "# Transverse "
 if( DEBUG ) print "# FreqRes_T_0, FreqRes_T = ", FreqRes_T_0, FreqRes_T 
 if( DEBUG ) print "# h_HOM_T, Mu_T, P_T = ", h_HOM_T, Mu_T, P_T ;
 if( DEBUG ) print "# R_T, Q_T = ", R_T, Q_T ; 

#***************************************************8
#
# Growth Rate of Longitudinal Coupled-Bunch Instabilities
#
#***************************************************8
if(DEBUG2 == 1){
 nres = int(FreqRes_L * Tref/NBunch) ;
 gratemax = 0.0
 for( mu=0; mu < NBunch;mu++)
 {
   ReZLxF[mu] = 0
   for(i=-nres*4;i<=nres*4;i++)
   {  
     f = 1/Tref * (i*NBunch + mu + NuS)
     ReZLxF[mu] = ReZLxF[mu] + f /( 1 + (Q_L*(FreqRes_L/f -f/FreqRes_L))^2 ) 
   }
   ReZLxF[mu] = R_L * ReZLxF[mu]
   GrowthRate_L[mu] = Current*Eta/(2*Eref_eV*NuS)*ReZLxF[mu]
   if( gratemax < GrowthRate_L[mu] )
   { 
       gratemax = GrowthRate_L[mu]
       mu_gratemax = mu 
   } 
 } 
 if( gratemax == 0 ){
      taumin = 0
 }else{
      taumin = 1./gratemax
 }
 print "# mu,Grate_L,Gtime_L= ", mu_gratemax, gratemax, taumin 
}
#***************************************************8
#
# Growth Rate of Transverse Coupled-Bunch Instabilities
#
#***************************************************8
if(DEBUG2 == 1){
 nres = int(FreqRes_T * Tref/NBunch) ;
 gratemax = 0.0
 for( mu=0; mu < NBunch;mu++)
 {
   ReZT[mu] = 0
   for(i=-nres*4;i<=nres*4;i++)
   {  
     f = 1/Tref * (i*NBunch - mu - NuBeta)
     ReZT[mu] = ReZT[mu] + (f/FreqRes_T) /( 1 + (Q_T*(FreqRes_T/f -f/FreqRes_T))^2 ) 
   }
   ReZT[mu] = R_T * ReZT[mu]
   GrowthRate_T[mu] = Current*Beta/(2*Eref_eV*Tref)*ReZT[mu]
   if( gratemax < GrowthRate_T[mu] )
   { 
       gratemax = GrowthRate_T[mu]
       mu_gratemax = mu 
   } 
 } 
 if( gratemax == 0 ){
      taumin = 0
 }else{
      taumin = 1./gratemax
 }
 print "# mu,Grate_T,Gtime_T= ", mu_gratemax, gratemax, taumin 
}
#***************************************************8
#
# ======== Get FreqResAcc 
#
#***************************************************8

RAcc0 = RperQAcc*QAcc0

I0 = Current 

Pb = I0 * ( U0 + Uid )
Pc = Vc0*Vc0 / (2.0*RAcc0)
BetaCOpt = 1 + Pb/Pc 

if(DEBUG) print "# RAcc0, QAcc0     = ", RAcc0, QAcc0
if(DEBUG) print "# Pb, Pc           = ", Pb, Pc ; 
if(DEBUG) print "# Beta Couple      = ", BetaC ;
if(DEBUG) print "# Beta Couple(Opt) = ", BetaCOpt ;

QAccL = QAcc0 / (1+BetaC) 
RAccL = RAcc0 / (1+BetaC) 

if(DEBUG) print "# RAccL, QAccL     = ", RAccL, QAccL

#***************************************************8
#
# get Psi  
#
#***************************************************8
cosPhi = U/Vc0
sinPhi = sqrt(1 - cosPhi*cosPhi)

Phi = atan2(sinPhi,cosPhi);

DPsi=-5 /360*2*PI

Theta = Phi ; 
TanPsi0 = - 2*I0*RAccL/Vc0 * sin(Theta) ;
Psi0 = atan2(TanPsi0,1) ;

Theta = Phi - DPsi 
TanPsi = 1/cos(DPsi)*( sin(DPsi) - 2*I0*RAccL/Vc0 * sin(Theta) ) ;
Psi = atan2(TanPsi,1.0) ;

if(DEBUG) print "# DPsi = Theta - Psi = ", DPsi/PI*180
if(DEBUG) print "# Phi, Theta = ", Phi/PI*180, Theta/PI*180 
if(DEBUG) print "# Psi0, Psi = ", Psi0/PI*180,Psi/PI*180 

Delta = -TanPsi/(2*QAccL)
FreqResAcc = FreqAcc/(1+Delta) 

Vb_re = - 2*I0*RAccL*cos(Psi)*cos(Psi)
Vb_im = - 2*I0*RAccL*cos(Psi)*sin(Psi)
Vb0 = sqrt( Vb_re*Vb_re + Vb_im*Vb_im )

Vc_re = Vc0 * cos(Phi)
Vc_im = Vc0 * sin(Phi)
Vc0 = sqrt( Vc_re*Vc_re + Vc_im*Vc_im )

Vg_re = Vc_re - Vb_re
Vg_im = Vc_im - Vb_im

Vg0 = sqrt( Vg_re*Vg_re + Vg_im*Vg_im )

if(DEBUG) print "# Vb0, Re[Vb], Im[Vb] = ", Vb0, Vb_re, Vb_im
if(DEBUG) print "# Vc0, Re[Vc], Im[Vc] = ", Vc0, Vc_re, Vc_im
if(DEBUG) print "# Vg0, Re[Vg], Im[Vg] = ", Vg0, Vg_re, Vg_im
if(DEBUG) print "# Delta, DFreq, FreqResAcc = ", Delta, Delta*FreqAcc, FreqResAcc 

#Pg = Vc0*Vc0*BetaC/(2*RAcc0)
#if(DEBUG) print "# Pg, Pc+Pb = ", Pg, Pb+Pc ; 
#Ig = 4*sqrt(Pg*BetaC/(2*RAcc0))
#Vg0 = Ig*RAcc0/(1+BetaC)

#***************************************************8
#
# ======== Get PhiSg, TauSg, NuSg ======
#            for Vc0 = Vg0 
# 
#    for Current = 0 (Initial State)
#
#***************************************************8

if(DEBUG) print "# Vg0 = ", Vg0 ;
sinPhiSg = U/Vg0
cosPhiSg = sqrt( 1 - sinPhiSg*sinPhiSg)
PhiSg = atan2(sinPhiSg,cosPhiSg)
OmegaSg = sqrt( Eta*Vg0/(Eref_eV*Tref)*OmegaAcc*cosPhiSg )
FreqSg = OmegaSg/(2*PI)
NuSg = FreqSg/Freqref
tauSg = PhiSg/OmegaAcc 
if(DEBUG) print "# I=0: NuSg, PhiSg(deg.) = ", NuSg, PhiSg/(2*PI)*360 


TauE = Eref_eV/U*Tref ;
if(DEBUG) print "# TauE = ", TauE ;
SettlingTime = TauE*2.0 ; 
RiseTime = TauE*2.0 ; 


Beam_nde     = 0.0 ; 
Beam_nde_rms = 0.0 ; 
Beam_tau     = tauSg 
Beam_tau_rms = Beam_nde_rms*Eta/OmegaS
Beam_a       = 0.0 
Beam_a_rms   = 1e-3
Beam_ap      = 0.0
Beam_ap_rms  = 1.952362e-05

Bunch_nde_rms = nde_rms_Natural ; 
Bunch_tau_rms = Bunch_nde_rms*Eta/OmegaS
Bunch_a_rms   = 1.0e-3 
Bunch_ap_rms  = 1.952362e-05
#***************************************************8
#
# Growth Rate of Coupled-Bunch Instabilities
#    by Acceleration Cavities
#
#***************************************************8

if(DEBUG2 == 1){
 print "# FreqResAcc = ", FreqResAcc
 print "# FreqResAcc - FreqAcc = ", FreqResAcc - FreqAcc
 print "# RAccL, QAccL     = ", RAccL, QAccL

 nres = int( FreqResAcc/(Freqref*NBunch) );
 for( mu=0 ; mu < NBunch ; mu++ )
 {
   ReZLxF[mu] = 0
   ImZLxF[mu] = 0
   for(i=-nres*2;i<=nres*2;i++)
   {  
     f = Freqref * (i*NBunch + mu + NuS)
     a =  f /( 1 + ( QAccL*( FreqResAcc/f -f/FreqResAcc ) )^2 ) 
     ReZLxF[mu] = ReZLxF[mu] + a 
     ImZLxF[mu] = ImZLxF[mu] + a *(-1)*QAccL*( FreqResAcc/f - f/FreqResAcc )
   }
   ReZLxF[mu] = RAccL * ReZLxF[mu]
   ImZLxF[mu] = RAccL * ImZLxF[mu]
   b = Current*Eta/(2*Eref_eV*NuS)
   GrowthRateAcc[mu] = b * ReZLxF[mu]
   OmegaShiftAcc[mu] = b * (-ImZLxF[mu])
   print "# mu,Grate_Acc,Fshift,Gtime_Acc=",mu,GrowthRateAcc[mu],OmegaShiftAcc[mu]/(2*PI),1/GrowthRateAcc[mu]
 } 
}

#***************************************************8
#
# ======== Get Instability Resonance Frequency
#
#***************************************************8

DeltaFreqRes=0

P=int(FreqRes0/Freqref/NBunch)
Mu=int(FreqRes0/Freqref)-P*NBunch
FreqRes=(P*NBunch+Mu+NuS)*Freqref + DeltaFreqRes

if(DEBUG) print "# P, Mu = ", P, Mu ;

CBISTOP=1

NumOfMonitorBunch=11
 MonBunchNo[0] = 5 
 MonBunchNo[1] = 10
 MonBunchNo[2] = 15
 MonBunchNo[3] = 20
 MonBunchNo[4] = 25
 MonBunchNo[5] = 30 
 MonBunchNo[6] = 35 
 MonBunchNo[7] = 40 
 MonBunchNo[8] = 45
 MonBunchNo[9] = 50
 MonBunchNo[10] = 55



randomseed = 101 

}
{
printf "%-15s%-15s%-15s%-15s\n", "NParticle", "NBunch", "NTurn", "SamplingStep"
printf "%-15d%-15d%-15d%-15d\n", NParticle, NBunch, NTurn, SamplingStep

printf "%-15s\n", "NBunchWithQ"
printf "%-15d\n", NBunchWithQ

printf "%-15s%-15s\n", "SettlingTime", "RiseTime"
printf "%-15.7e%-15.7e\n", SettlingTime, RiseTime 

printf "%-15s%-15s\n", "Eref_eV", "Tref" 
printf "%-15.7e%-15.7e\n",Eref_eV,  Tref

printf "%-15s%-15s%-15s\n", "NuBeta", "Eta", "Je"
printf "%-15.7e%-15.7e%-15.7e\n", NuBeta, Eta, Je 

printf "%-15s%-15s%-15s%-15s\n", "Chrom0", "Chrom1", "ChromFreq", "ChromPhase"
printf "%-15.7e%-15.7e%-15.7e%-15.7e\n", Chrom0, Chorm1, ChromFreq, ChromPhase  

printf "%-15s%-15s\n", "U0(eV)", "Uid(eV)"
printf "%-15.7e%-15.7e\n", U0, Uid

printf "%-15s%-15s\n", "nde_rms_Natural", "Emit_Natural"
printf "%-15.7e%-15.7e\n", nde_rms_Natural, Emittance_Natural  

printf "%-15s%-15s\n", "LDamping", "LExcitation"
printf "%-15d%-15d\n", LDamping, LExcitation

printf "%-15s%-15s\n", "TDamping", "TExcitataion"
printf "%-15d%-15d\n", TDamping, TExcitation 

printf "%-15s\n", "Current"
printf "%-15.7e\n", Current

printf "%-15s%-15s%-15s%-15s\n", "Beam_nde", "Beam_nde_rms", "Beam_tau", "Beam_tau_rms"
printf "%-15.7e%-15.7e%-15.7e%-15.7e\n", Beam_nde, Beam_nde_rms, Beam_tau, Beam_tau_rms ; 

printf "%-15s%-15s%-15s%-15s\n", "Beam_a", "Beam_a_rms", "Beam_ap", "Beam_ap_rms"
printf "%-15.7e%-15.7e%-15.7e%-15.7e\n", Beam_a, Beam_a_rms, Beam_ap, Beam_ap_rms ; 

printf "%-15s%-15s\n", "Bunch_nde_rms", "Bunch_tau_rms"
printf "%-15.7e%-15.7e\n", Bunch_nde_rms, Bunch_tau_rms ; 

printf "%-15s%-15s\n", "Bunch_a_rms", "Bunch_ap_rms"
printf "%-15.7e%-15.7e\n", Bunch_a_rms, Bunch_ap_rms ; 

printf "%-15s\n", "NAcc"
printf "%-15d\n", NAcc

printf "%-15s%-15s%-15s\n", "RperQ", "Ql", "FreqRes"
printf "%-15.7e%-15.7e%-15.7e\n", RperQAcc, QAccL, FreqResAcc

printf "%-15s%-15s%-15s\n", "Vg0","Harmonics", "PhaseAcc" ;
printf "%-15.7e%-15d%-15.7e\n", Vg0, Harmonics, PhaseAcc ; 

printf "%-15s\n", "NLRes"
printf "%-15d\n", NLRes

if( NLRes > 0 ){
  printf "%-15s%-15s%-15s\n", "RperQ_L", "Q_L", "FreqRes_L"
  printf "%-15.7e%-15.7e%-15.7e\n", RperQ_L, Q_L, FreqRes_L
}

printf "%-15s\n", "NTRes"
printf "%-15d\n", NTRes

if( NTRes > 0 ){
  printf "%-15s%-15s%-15s%-15s%-15s\n", "RperQ_T", "Q_T", "FreqRes_T", "Beta", "Alpha"
  printf "%-15.7e%-15.7e%-15.7e%-15.7e%-15.7e\n", RperQ_T, Q_T, FreqRes_T, Beta, Alpha
}

printf "%-15s\n", "NLDevice"
printf "%-15d\n", NLDevice

printf "%-15s\n", "NTDevice"
printf "%-15d\n", NTDevice

printf "%-15s\n", "CBISTOP"
printf "%-15d\n", CBISTOP 

printf "%-15s\n", "NumOfMoniorBunch" 
printf "%-15d\n", NumOfMonitorBunch


if( NumOfMonitorBunch > 0 )
{
  printf "%s\n", "NumOfMoniorBunchlist" 
  for(i=0; i<NumOfMonitorBunch;i++)
  {
    printf "%-5d ", MonBunchNo[i] 
  }
  printf "\n" ;
}


printf "%-15s\n", "randomseed"
printf "%-15d\n", randomseed
}' > $outfile