#!/bin/bash
############################################################################
#script to submit QCT calculations using the QCTtraj code
#
#	QCTraj uses as input V0 [m/s]  alpha [m/s]  Tn[k]  D2/H2 

#arrays to hold beam parameters
v0=();		alpha=();		Tn=();		E0=();		Ei=();		DorH=(); 	S0exp=();		

#Seeded Molecular H2 beams Tsurf =120K --> Physical chemistry chemical physics PCCP, 12, 6499, (2010) 				Aray #	Exp data origin
Tn+=("1740");		v0+=("3923");		alpha+=("1105");		DorH+=("H");		S0exp+=("0.00145");	#0 	alm exp Cu111
Tn+=("1740");           v0+=("4892");           alpha+=("1105");                DorH+=("H");            S0exp+=("0.00400");	#1 	alm exp Cu111
Tn+=("1740");           v0+=("5906");           alpha+=("945");                 DorH+=("H");            S0exp+=("0.01010");	#2 	alm exp Cu111
Tn+=("2000");           v0+=("3857");           alpha+=("995");                 DorH+=("H");            S0exp+=("0.00283");	#3 	alm exp Cu111
Tn+=("2000");           v0+=("4625");           alpha+=("1032");                DorH+=("H");            S0exp+=("0.00556");	#4 	alm exp Cu111
Tn+=("2000");           v0+=("6432");           alpha+=("886");                 DorH+=("H");            S0exp+=("0.02606");	#5 	alm exp Cu111

#Pure H2 Molecular beams Tsurf =120K --> Physical chemistry chemical physics PCCP, 12, 6499, (2010)
Tn+=("1435");           v0+=("5417");           alpha+=("826");                 DorH+=("H");            S0exp+=("0.00197");     #6	alm exp Cu111
Tn+=("1465");           v0+=("5446");           alpha+=("830");                 DorH+=("H");            S0exp+=("0.00197");     #7	alm exp Cu111
Tn+=("1740");           v0+=("5906");           alpha+=("945");                 DorH+=("H");            S0exp+=("0.01023");     #8	alm exp Cu111
Tn+=("1855");           v0+=("6139");           alpha+=("899");                 DorH+=("H");            S0exp+=("0.01714");     #9	alm exp Cu111
Tn+=("2000");           v0+=("6431");           alpha+=("886");                 DorH+=("H");            S0exp+=("0.02624");     #10	alm exp Cu111
Tn+=("2100");           v0+=("6674");           alpha+=("913");                 DorH+=("H");            S0exp+=("0.04918");     #11	alm exp Cu111
Tn+=("2300");           v0+=("6590");           alpha+=("1351");                DorH+=("H");            S0exp+=("0.05995");     #12	alm exp Cu111

#Seeded D2 Molecular beams Tsurf =120K --> Physical chemistry chemical physics PCCP, 12, 6499, (2010)
Tn+=("2100");           v0+=("3925");           alpha+=("816");                 DorH+=("D");            S0exp+=("0.02109");     #13      alm exp Cu111
Tn+=("2100");           v0+=("4595");           alpha+=("782");                 DorH+=("D");            S0exp+=("0.03854");     #14      alm exp Cu111
Tn+=("2100");           v0+=("5377");           alpha+=("829");                 DorH+=("D");            S0exp+=("0.14036");     #15      alm exp Cu111
Tn+=("2100");           v0+=("5658");           alpha+=("866");                 DorH+=("D");            S0exp+=("0.25818");     #16      alm exp Cu111
Tn+=("2100");           v0+=("6132");           alpha+=("849");                 DorH+=("D");            S0exp+=("0.33818");     #17      alm exp Cu111

#Pure D2 Molecular beams Tsurf =120K -->Physical chemistry chemical physics PCCP, 12, 6499, (2010)
Tn+=("1435");           v0+=("4014");           alpha+=("299");                 DorH+=("D");            S0exp+=("0.00177");     #18      alm exp Cu111
Tn+=("1790");           v0+=("4196");           alpha+=("614");                 DorH+=("D");            S0exp+=("0.00884");     #19      alm exp Cu111
Tn+=("1670");           v0+=("4337");           alpha+=("371");                 DorH+=("D");            S0exp+=("0.00733");     #20      alm exp Cu111
Tn+=("1905");           v0+=("4374");           alpha+=("685");                 DorH+=("D");            S0exp+=("0.01371");     #21      alm exp Cu111
Tn+=("1975");           v0+=("4461");           alpha+=("687");                 DorH+=("D");            S0exp+=("0.02305");     #22      alm exp Cu111

#H2 + Cu111 --> Kun from Ludo's group experiments <UNPUBLISHED>
Tn+=("1173.15");        v0+=("4714.489");       alpha+=("1526.307");            DorH+=("H");            S0exp+=("0.01399");     #23	 Ludo Cu111
Tn+=("1223.15");        v0+=("4853.248");       alpha+=("1559.036");            DorH+=("H");            S0exp+=("0.01814");     #24	 Ludo Cu111
Tn+=("1273.15");        v0+=("4553.746");       alpha+=("1834.320");            DorH+=("H");            S0exp+=("0.02313");     #25	 Ludo Cu111
Tn+=("1323.15");        v0+=("4800.152");       alpha+=("1771.183");            DorH+=("H");            S0exp+=("0.01150");     #26	 Ludo Cu111

#D2 + Cu111 --> Kun from Ludo's group experiments <UNPUBLISHED>
Tn+=("1373.15");        v0+=("2988.02");        alpha+=("998.7326");            DorH+=("D");            S0exp+=("0.02121");     #27	 Ludo Cu111
Tn+=("1273.15");        v0+=("3091.629");       alpha+=("960.4899");            DorH+=("D");            S0exp+=("0.01749");     #28	 Ludo Cu111
Tn+=("1173.15");        v0+=("3133.507");       alpha+=("881.5357");            DorH+=("D");            S0exp+=("0.01340");     #29	 Ludo Cu111
Tn+=("1473.15");        v0+=("3256.816");       alpha+=("1133.679");            DorH+=("D");            S0exp+=("0.02332");     #30	 Ludo Cu111

#D2 + Cu111 --> From my 3th author paper with gernot and kun
Tn+=("1173.15");        v0+=("3216.59764");     alpha+=("1155.62446");          DorH+=("D");            S0exp+=("0.01287");     #31      Ludo Cu111
Tn+=("1223.15");        v0+=("3408.985424");    alpha+=("1150.780309");         DorH+=("D");            S0exp+=("0.01485");     #32      Ludo Cu111
Tn+=("1223.15");        v0+=("3496.892782");    alpha+=("1227.993863");         DorH+=("D");            S0exp+=("0.02073");     #33      Ludo Cu111
Tn+=("1273.15");        v0+=("3518.912431");    alpha+=("1463.893456");         DorH+=("D");            S0exp+=("0.02430");     #34      Ludo Cu111

#D2 + Cu211 --> From my 3th author paper with gernot and kun
Tn+=("1223.15");        v0+=("3172.972");       alpha+=("1136.581");            DorH+=("D");            S0exp+=("0.00738");     #35      Ludo Cu111
Tn+=("1173.15");        v0+=("3211.937");       alpha+=("1281.84");             DorH+=("D");            S0exp+=("0.00829");     #36      Ludo Cu111
Tn+=("1373.15");        v0+=("3220.79");        alpha+=("1092.302");            DorH+=("D");            S0exp+=("0.00976");     #37      Ludo Cu111
Tn+=("1273.15");        v0+=("3314.637");       alpha+=("1315.568");            DorH+=("D");            S0exp+=("0.00916");     #38      Ludo Cu111

#Table S6 Diaz 2009 SRP Cu111 SOM, Berger /Rendulic broad beams --> <Ei> = 2.7K_B T_N
Tn+=("1118.07");        v0+=("3501");           alpha+=("1996");                DorH+=("H");            S0exp+=("0.00916");     #39      
Tn+=("1331.89");        v0+=("3556");           alpha+=("2342");                DorH+=("H");            S0exp+=("0.00916");     #40
Tn+=("1438.83");        v0+=("3381");           alpha+=("2611");                DorH+=("H");            S0exp+=("0.00916");     #41
Tn+=("1501.19");        v0+=("3152");           alpha+=("2819");                DorH+=("H");            S0exp+=("0.00916");     #42
Tn+=("1581.35");        v0+=("3219");           alpha+=("2903");                DorH+=("H");            S0exp+=("0.00916");     #43


#################################
#constansts
massH=1.0072 
massD=2.0141

Eshift="-0.0"


#for loop over the molecular beam parameters
for i in {0..43}
do


  #Calculate E0 [eV] from v0 [m/s]
  if [ ${DorH[$i]} == "H" ]; then
  
    E0[$i]=`echo "scale=5; 0.5 * 2 * $massH * ${v0[$i]} * ${v0[$i]} * 0.000000010364269 / 1" | bc -l`
    #subtract Eshift from E0
    E0[$i]=`echo " ${E0[$i]} ${Eshift} " | bc -l`
    #recalculate v0 based on new E0
    v0[$i]=`echo "scale=5; sqrt( (2 * ${E0[$i]} ) / ( 2 * $massH ) ) * 9822.695 / 1" | bc -l`
  elif [ ${DorH[$i]} == "D" ]; then 
    E0[$i]=`echo "scale=5; 0.5 * 2 * $massD * ${v0[$i]} * ${v0[$i]} * 0.000000010364269 / 1" | bc -l`
    E0[$i]=`echo " ${E0[$i]} ${Eshift} " | bc -l`
    v0[$i]=`echo "scale=5; sqrt( (2 * ${E0[$i]} ) / ( 2 * $massD ) ) * 9822.695 / 1" | bc -l`
  fi
  
  echo ${Tn[$i]} ${v0[$i]} ${E0[$i]} ${alpha[$i]} ${Ei[$i]}
done

mkdir BEAMS
echo "# table with data for H2 + Cu  " > table.dat
echo "# T_N [K] | V0 [m/s] | E0 [m/s] | alpha [m/s] | S0 exp | S0 calc " >> table.dat 
for i in 0 1 2 3 4 5 6 7 8 9 10 11 12 23 24 25 26 39 40 41 42 43
do
./analysis ${Tn[$i]} ${E0[$i]} ${alpha[$i]} 0.001 
cp state.occupencies BEAMS/beam${i}.state.occupencies
cp out.info BEAMS/beam${i}.info
S0calc=`tail -1 out.info | awk '{print $3}'`
rm out.info state.occupencies
echo "${Tn[$i]} ${v0[$i]} ${E0[$i]} ${alpha[$i]} ${S0exp[$i]} $S0calc" >> table.dat
done #i
cat table.dat