"GENERAL HO;

f[n_,x_]:= Abs[((1/Pi)^(1/4)HermiteH[n,x])/(E^(x^2/2)Sqrt[2^n n!])]^2

Plot[Evaluate@
  Append[Table[f[n,x]+n+1/2,{n,0,7}],x^2/2],{x,-4,4},Filling->Table[n->n-1/2,{n,1,8}]]

" this is not real world picture,
as long as r correctnot realistic if not r correct";

" If wave packets than overlaps or squared functions for 
transitions and steady states, rspktvl'

"MORSE;

"Boltzman * FC ;

"see for I2 - you can do for HCl co2 etc.;

"two states B and X ground, see hertzberg;

"units  r-angstrem, w cm-1,D cm-1 mas kg;c cm

Clear[x,b]

r[x]=2.6663
r[b]=3.0247
w[x]=214.502
w[b]=125.697
de[x]=12439
de[b]=5168.0

mu=0.1269/2

c=2.9979 10^10

h=6.626 10^-34
k=1.381 10^-23
Na=6.022 10^23
Pi
cont=  2 Pi (Sqrt[((c mu)/(h Na))])10^-10

re=2.6663 10^-10

" this is all in SI;

νe=0.2145 10^5

masp=1836 *9.1 10^-31

μ=(1/((1/(126.9 *masp))+(1/(126.9*masp))))

fre = νe /100

Solve[fre==(Sqrt[(kk/μ)])*(8065.5*h/(1.6 10^-19))/(2 Pi),kk]

(de[x]/8065.5)

dd=(de[x]/8065.5)*1.6 (10^-19)

νe=(Sqrt[(k/μ)])

β=2 Pi (c/100) νe  (Sqrt[(μ/(2  dd))])

Plot[dd  (1-  Exp[-β  (r-re)])^2,{r,2.1 (10^(-10)),4.1 (10^(-10))}]

NOW EXXCITED electronic state

r[b]

reb=3.0247 10^-10

" this is all in SI;

νeb=0.125 10^5

masp=1836 *9.1 10^-31

μ=(1/((1/(126.9 *masp))+(1/(126.9*masp))))

freb = νeb  /100

Solve[freb==(Sqrt[(kkb/μ)])*(8065.5*h/(1.6 10^-19))/(2 Pi),kkb]

(de[b]/8065.5)

ddb=(de[b]/8065.5)*1.6 (10^-19)

βb=2 Pi (c/100) νeb  (Sqrt[(μ/(2  ddb))])

Plot[ddb  (1-  Exp[-βb (r-reb)])^2,{r,2.5 (10^(-10)),4.1 (10^(-10))}]

Plot[{(5*10^-19)+ddb  (1-  Exp[-βb (r-reb)])^2,dd  (1-  Exp[-β  (r-re)])^2},{r,2.1 (10^(-10)),4.1 (10^(-10))}]

d[x]=2 de[x] /w[x]

d[b]=2 de[b]/w[b]

beta[x]=cont r[x] (Sqrt[w[x]])
beta[b]=cont r[b](Sqrt[w[b]])

alfa[x]=beta[x]/(Sqrt[d[x]])

alfa[b]=beta[b]/(Sqrt[d[b]])

"HO;

H[t_,e_,v_]:=((HermiteH[v,z=beta[e] ((t/(r[e]))-1)] Exp[- (z^2)/2])/(Sqrt[(((Sqrt[Pi ])  (2^v)  v! r[e])/(beta[e]))]))

"see paper jpc 95 (21) 8019;

vx=0
vb=34

plot[f_]:= Plot[{f[t,x,vx]^2, 10+f[t,b,vb]^2,22+f[t,x,vx]*f[t,b,vb]^2},{t,2.4,3.8},PlotRange->{0,30},PlotPoints->100,Frame->True,   FrameTicks->{Automatic,None},FrameLabel->{ "t (Angstroms)"," X*X B*B X*B  "},DisplayFunction->Identity]

Show[GraphicsArray[{plot[H]}],DisplayFunction->$DisplayFunction]

alfa[x]

"MORSE has its own solutions eigenfunctions";

M[t_,e_,v_]:= (LaguerreL[v,a=2 d[e]-2 v -1 ,y=2  d[e]  Exp[-alfa[e] ((t/(r[e]))-1)]] (Sqrt[Exp[-y] (y^a)]))/((Sqrt[((((∑_(s=0)^v)((Gamma[s+a])/(s!)) ) r[e])/(alfa[e]))]))

"not;

M[t_,e_,v_]:= (LaguerreL[v,a=2 d[e]-2 v -1 ,y=2  d[e]  Exp[-alfa[e] ((t/(r[e]))-1)]] (Sqrt[Exp[-y]  (y^a)]))/((Sqrt[((((∑_(s=0)^v)(((Gamma[s+a])/(s!))) ) r[e])/(alfa[e]))]))

plot[f_]:= Plot[{f[t,x,vx]^2, 10+f[t,b,vb]^2,22+f[t,x,vx]*f[t,b,vb]},{t,2.4,4.8},PlotRange->{0,30},PlotPoints->100,Frame->True,   FrameTicks->{Automatic,None},FrameLabel->{ "t (Angstroms)"," X*X B*B X*B  "},DisplayFunction->Identity]

Show[GraphicsArray[{plot[H],plot[M]}],DisplayFunction->$DisplayFunction]

"Now take care of Botznman;

T=298

w[x]

BF[v_]:= Exp[-v w[x] h c /(k T)]

BF[0]

BF[34]

NIntegrate[M[t,x,vx] M[t,b,vb], {t,2.4,2.9}]^2

"FCF[vx_,vb_,f_]:=NIntegrate[f[t,x,vx] f[t,b,vb], {t,2.4,2.9},MaxRecursion->20, AccuracyGoal->8]^2

FCF[vx_,vb_,f_]:=NIntegrate[f[t,x,vx] f[t,b,vb], {t,2.4,2.9},MaxRecursion->20, AccuracyGoal->8]^2

RI[vx_,vb_,f_]:=BF[vx]  FCF[vx,vb,f]

vx

plot[f_]:=Show[Table[ListPlot[Table[{vb,RI[vx,vb,f]},{vb,0,34}],PlotJoined->True,PlotRange->All,DisplayFunction->Identity],{vx,0,2}],Frame->True,FrameLabel->{" v for B state", "relative intensity RI"}]

T

vx

k

Mx0[t_]:= (LaguerreL[0,a=2 d[x]-2 0-1 ,y=2  d[x]  Exp[-alfa[x] ((t/(r[x]))-1)]] (Sqrt[Exp[-y] (y^a)]))/((Sqrt[((((∑_(s=0)^0)((Gamma[s+a])/(s!)) ) r[x])/(alfa[x]))]))

Mx0[2.]

Mb34[t_]:= (LaguerreL[34,a=2 d[b]-2 34 -1 ,y=2  d[b]  Exp[-alfa[b] ((t/(r[b]))-1)]] (Sqrt[Exp[-y] (y^a)]))/((Sqrt[((((∑_(s=0)^34)((Gamma[s+a])/(s!)) ) r[b])/(alfa[b]))]))

Mb34[2.]

ht=Table[Mx0[t] Mb34[t], {t,2.4,2.9,0.002}];

"for vx=0

Exp[-vx w[x] h c /(k T)]

"for vx=1

Exp[-2 w[x] h c /(k T)]

allt=Table[ 1  Mx0[t] Mb34[t]  Mx0[t] Mb34[t], {t,2.4,2.9,0.002}];

ListPlot[{ht,allt},PlotJoined->True,PlotRange->All]

ListPlot[ht*ht*,PlotJoined->True,PlotRange->All]

Show[GraphicsArray[{plot[H],plot[M]}],PlotRange->{0,0.2},DisplayFunction->$DisplayFunction]

Show[GraphicsArray[{plot[M]}],DisplayFunction->$DisplayFunction]

(1./5210)*10^8

(1./11600)*10^8

(1./7900)*10^8

19193.8579654510-8620.6896551724

de[x]

de[b]

T=298

"so we know the target wave functions, what is the EM field that would drive system in those selected states (quantum control or selection, event selection0;

" use an electric field ;

EMF[t_]:= CC Exp[- aa (t-t1)^2]  Cos[ww t +fi]

CC=1

aa=1

t1=2

ww=4

fi=Pi

Plot[EMF[t],{t,-10,10},PlotRange->All]

"-1 to 5"

6./250

tempEMF=Table[EMF[i],{i,-1,5,0.025}]

FourierTransform[1/(t^2+1),t,ω]

InverseFourierTransform[%, ω,t]

FourierTransform[1,t,ω]

FourierTransform[ CC Exp[- aa (t-t1)^2]  Cos[ww t +fi],t,ω]

frEMF=Abs[Fourier[tempEMF]]

q2=Chop[%]

ListPlot[q2]