{VERSION 3 0 "IBM INTEL NT" "3.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 }{CSTYLE "" -1 256 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 258 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 261 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 } {CSTYLE "" -1 264 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Title" 0 18 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 1 0 0 0 0 0 0 }3 0 0 -1 12 12 0 0 0 0 0 0 19 0 }{PSTYLE "Author" 0 19 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 8 8 0 0 0 0 0 0 -1 0 }} {SECT 0 {EXCHG {PARA 18 "" 0 "" {TEXT -1 94 "Symbolic Programmes for A nalysis of Vibrational and Rotational Spectra of Diatomic Molecules \+ " }}{PARA 0 "" 0 "" {TEXT -1 79 " This work sheet in Maple con tains procedures in four suites described " }}{PARA 19 "" 0 "" {TEXT -1 3 "by " }{TEXT 263 33 "F. M. Fernandez and J. F. Ogilvie" }}{PARA 0 "" 0 "" {TEXT -1 11 "in journal " }{TEXT 264 9 "MapleTech" }{TEXT -1 390 ", volume 5, 1998, No. 1, pages 42 - 46. These procedures for \+ symbolic computation in Maple (version V, release 5.1) allow one to fo rm the principal expressions required to undertake a comprehensive ana lysis of vibration-rotational spectra of diatomic molecules in an elec tronic state singlet-SIGMA-+ according to conventional analytic treatm ent in the spirit of Dunham (cf. J. F. Ogilvie, " }{TEXT 256 65 "The V ibrational and Rotational Spectrometry of Diatomic Molecules" }{TEXT -1 1018 ", Academic Press, London, 1998). To relate to measurements \+ of frequencies or wavenumbers of spectral transitions, two procedures \+ form, in the first suite, expressions Y_kl for spectral term coefficie nts associated with mechanical effects, i.e. vibrational and rotationa l motions of atomic centres about the centre of molecular mass, and, i n the second suite, expressions Z_kl for further spectral term coeffic ients associated with extra-mechanical effects, i.e. in relation to fa ilure of electrons to follow perfectly one or other nucleus. To relat e to measurements of intensities of spectral lines, two additional pro cedures form, in the third suite, vibrational matrix elements and, in the fourth suite, the Herman-Wallis factor F_0^v' that take s into account vibration-rotational interactions that affect relative \+ intensities within bands. Use of these suites is described specifical ly in the paper in MapleTech mentioned above; their applications are d escribed at length in the book specified above." }}}{EXCHG {PARA 0 "" 0 "" {TEXT 265 50 "1 Term coefficients Y_kl for mechanical effects " }{TEXT -1 2 " " }}{PARA 0 "" 0 "" {TEXT -1 342 " The first sui te DUNHAM of procedures below generates Dunham's term coefficients Y[k ,l] for vibration-rotational energies of diatomic molecules in electro nic state singlet-SIGMA+ with the method HPM (TP-hypervirial-Hellmann- Feynman perturbation theory); the quantity JB is considered to be inde pendent of perturbation parameter lambda; " }}{PARA 0 "" 0 "" {TEXT -1 2028 " definitions: \+ lambda^2=2*Be/We=gamma=g, vh=v+1/2, JB=JJ*g/2, JJ=J (J+1).\n To use procedure tp1 to generate term coefficients Y[k,l ,m], in which k is the power of (v+1/2), l is the power of [J(J+1)] an d m is the order of the contribution (m=1 for the leading contribution to Y[k,l], m=2 for the first correction etc.), invoke \"tp1(p1):\" in which the argument \"p1\" is the required order of perturbation theor y. This procedure should work with values of p1 up to at least p1=14. All expressions Y[k,l,m] according to the selected order of perturba tion theory are automatically generated, and are invoked simply on typ ing \"Y[k,l,m];\" with appropriate values of k, l and m. For p1=14, t he range of Y[k,l,1] is from k=8,l=0 with a_1^14 ... a_14 to k=0,l=8 w ith a_1^6 ... a_6.\n To use procedure tp2 to generate term coeffi cients Y[k,l,m], with the same meaning of symbols as for tp1, invoke \+ \"tp2(p2):\" in which the argument \"p2\" is the required order of per turbation theory. This procedure should work with values of p2 up to a t least p2=20, but the duration of execution of tp2 for a given order \+ is about three times that of tp1. For p2=20, the range of Y[k,l,1] is from k=11,l=0 with a_1^20 ... a_20 to k=0,l=11 with a_1^9 ... a_9.\n \+ By default expressions Y[k,l,m] appear in terms of primary rotation al parameter Be, reciprocal W1 of primary vibrational parameter (harmo nic vibrational wavenumber) We, and coefficients a_j in Dunham's funct ion V(x) for potential energy with x=(R-Re)/Re. If these expressions \+ are preferred in terms of coefficients b_j of V(y) with y=(R-Re)/R or \+ c_j of V(z) with z=2(R-Re)/(R+Re) instead of a_j, then after the expre ssions Y[k,l,m] are formed with tp1 or tp2 call procedure \"converta(l ):\" with argument l as either b or c as required. To ensure that res ulting expressions Y[k,l,m] in terms of b_j or c_j appear in an attrac tive form, type an instruction \"sort(simplify(Y[k,l,m]));\" to displa y the desired expression." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 813 "restart: tp1:=proc(p1) # perturbation theory of order p1 \nlocal i,s,N,Q,p: global En,JJ,vh:\nif type(p1,odd) then p:=p1-1 else p:=p1 fi:\nQ[0,0]:=1: En[0]:=vh:\nfor s from 0 to p do\n if s> 0 the n\n Q[0,s]:=0:\n if type(s,odd) then En[s]:=0 else\n En[s]: =expand(1/s*add(i*(a.i/2*Q[i+2,s-i]\n +JJ*(-1)^i *(i+1)*Q[i,s-i]),i=1..s)):\n fi:\n fi:\n for N from 0 to p-s+1 do \n if type(N+1+s,odd) then Q[N+1,s]:=0 else\n Q[N+1,s]:=expand (1/(N+1)*(trunca(N-3)*N*(N-1)*(N-2)*Q[N-3,s]/4\n +2*N*E n[0]*trunca(N-1)*Q[N-1,s]\n +add(2*N*En[i]*trunca(N-1)* Q[N-1,s-i]\n -(N+i/2+1)*a.i*Q[N+i+1,s-i]\n \+ -(-1)^i*(i+1)*(2*N+i)*JJ*Q[N+i-1,s-i],i=1..s))):\n fi:\n od:\no d:\nEn[0]:=En[0]+JJ;\nfor i from 0 to p by 2 do xtrct(i) od; gc() end :" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "trunca:=proc(n) if n<0 then 0 else 1 fi: end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 264 "xtrct:=proc(n) local k,l,p; global En,Y,JJ,vh,Be,W1;\nfor k from 0 to n/2+1 do p:=subs(vh=0,En[n]); En[n]:=expand((En[n]-p)/vh);\n for l f rom 0 to n/2+1 do Y[k,l,3/2+n/4-k/2-l/2]:=combine(\n (Be*W1)^(l+n/2 )/W1*sort(2^(n/2)*coeff(p,JJ,l))) od\n od; En[n]:=0\nend:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 299 "converta:=proc(l) local j,k ; global b,c;\nif l=c then\n for k from 1 to 22 do a.k:=(-1)^k*(k+1)/ 2^k+mysum('(-1)^j*(k+1)!*\n c.(k-j)/(2^j*j!*(k+1-j)!)', 'j'=0..k- 1) od\nelif l=b then\n for k from 1 to 22 do a.k:=(-1)^k*(k+1)+mysum( '(-1)^j*(k+1)!*\n b.(k-j)/(j!*(k+1-j)!)','j'=0..k-1) od\nfi end: " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1407 "tp2:=proc(p2) \+ # perturbation theory of order p2\nlocal i,k,l,k1,l1,s,N,p: glob al e, Q:\nif type(p2,odd) then p:=p2-1 else p:=p2 fi:\nfor s from 0 to p do\n for N from 0 to p-s+1 do\n for k from 0 to p do\n for l from 0 to p do\n Q[N,k,l,s]:=0: e[k,l,s]:=0:\n od:\n \+ od:\n od:\nod:\nQ[0,0,0,0]:=1: e[0,0,0]:=1:\nfor s from 0 to p do\n \+ for k from 0 to s/2+1 do\n for l from 0 to s/2+1 do\n if s> 0 then\n if type(s,odd) then e[k,l,s]:=0 else\n e[k,l,s]:=s implify(1/s*add(i*(a.i/2*trunca(s/2+1-k-l)\n *Q[i+2,k ,l,s-i]\n +(-1)^i*(i+1)*trunca(l-1)*\n \+ trunca(s/2-k-l+1)*Q[i,k,l-1,s-i]),i=1..s)):\n fi:\n \+ fi:\n od:\n od:\n for N from 0 to p-s+1 do\n for k fr om 0 to (N+1+s)/2 do\n for l from 0 to (N+1+s)/2 do\n \+ if type(N+1+s,odd) then Q[N+1,k,l,s]:=0 else\n Q[N+1,k, l,s]:=simplify(1/(N+1)*(trunca(N-3)*N*(N-1)*(N-2)\n *Q[N- 3,k,l,s]/4\n +2*N*trunca(N-1)*trunca(k-1)*Q[N-1,k-1,l,s] \n +add(2*N*add(add(e[k1,l1,i]*\n \+ trunca(N-1)*Q[N-1,k-k1,l-l1,s-i],\n l1=0..l), k1=0..k)\n -(N+i/2+1)*a.i*Q[N+i+1,k,l,s-i]\n \+ -(-1)^i*(i+1)*(2*N+i)*trunca(l-1)*Q[N+i-1,k,l-1,s-i],i=1..s))):\n \+ fi\n od\n od\n od\nod: e[0,1,0]:=1: e[1, 0,0]:=1: extrct(p); gc()\nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 233 "extrct:=proc(n) local k,l,m; global e,Y,Be,W1;\nfor k from 0 \+ to n/2+1 do\n for l from 0 to n/2+1 do\n for m from 0 to n do\n \+ Y[k,l,3/2+m/4-k/2-l/2]:=combine((Be*W1)^(l+m/2)/W1*\n sort(2 ^(m/2)*e[k,l,m]))\n od\n od\nod\nend:" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 36 "tp1(12): Y[2,1,1]; # for instance" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 258 57 "2 Term coefficients Z_kl for extra-m echanical effects" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 915 " \+ Procedures in this suite ADNONAD allow one to generate additional t erm coefficients Z[k,l] to take account of adiabatic and nonadiabatic \+ contributions to vibration-rotational energies of diatomic molecules i n electronic state singlet-SIGMA+\n To use procedure tp3 to gener ate term coefficients Z[k,l], in which k is the power of (v+1/2) and l is the power of [J(J+1)], invoke \"tp3(P3):\", in which the argument \+ P3 is the required order of perturbation theory. This procedure shoul d work with values of P1 up to at least 12. All expressions Z[k,l] ar e automatically generated, and are displayed simply on typing \"Z[k,l] ;\" with appropriate values of k and l. For P1=12, the range of Z[k,l ] is k+l<=6.\n To convert coefficients a_j, K_j, q_j and d_j in V( x) etc. to c_j, u_j, t_j and s_j in V(z) etc. invoke \"convertacstu(): \" after running tp3; then display an expression with \"sort(simplify( Z[k,l]));\"." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1631 "restart: \+ tp3:=proc(P3) # perturbation theory of order P3\nlocal l,n, p,i,k,p1,la,P,s1,s2,s3,JB,ad1,ad2,ad3;\nglobal lista,Be,DE,En,JJ,We,en ,g,Q,vh;\nif type(P3,odd) then P:=P3-1 else P:=P3 fi;\nlista:=[g,Be,K. (1..12),q.(0..8),d.(0..8),a.(1..14)];\nQ[0,0]:=1; en[0]:=vh;\nfor p fr om 0 to P do\nif p>0 then\n Q[0,p]:=0:\n if type(p,odd) then en[p]:= 0 else\n en[p]:=expand(1/p*add(i*(a.i/2*Q[i+2,p-i]\n +JB*(-1) ^i*(i+1)*Q[i,p-i]),i=1..p)):\n fi:\nfi;\nfor n from 0 to P-p+1 do\n \+ if type(n+1+p,odd) then Q[n+1,p]:=0 else\n Q[n+1,p]:=expand(1/( n+1)*(trunca(n-3)*n*(n-1)*(n-2)*Q[n-3,p]/4\n +2*n*en[0] *trunca(n-1)*Q[n-1,p]\n +add(2*n*en[i]*trunca(n-1)*Q[n- 1,p-i]\n -(n+i/2+1)*a.i*Q[n+i+1,p-i]\n - (-1)^i*(i+1)*(2*n+i)*JB*Q[n+i-1,p-i],i=1..p))):\n fi: od; od; JB:=J J/2*g:\n# Form corrections adiabatic and nonadiabatic. \nla:=g^(1/2); \+ s1:=0; s2:=0; s3:=0;\nfor k from 0 to P do\n# correction nonadiabatic \+ vibrational\n ad1[k]:=expand(k*(k-1)/4*add(Q[k-2,p1]*la^(p1+k),p1=0.. P-k)\n +1/2/(k+1)*(1/2*add((i+2)*a.i*add(Q[i+k+2,p1]*la^(p1+k+i),\n p1=0..P-k-i),i=1..P-k)+JB*add((-1)^i*i*(i+1)\n *add(Q[i+k,p1]*l a^(p1+k+i),p1=0..P-k-i),\n i=1..P-k)+add(Q[k+2,p1]*la^(k+p1),p1=0.. P-k)));\n# correction nonadiabatic rotational\n ad2[k]:=expand(JB*add ((-1)^i*(i+1)*la^(i+k)\n *add(Q[i+k,l]*la^l,l=0..P-k-i),i=0..P-k ));\n# correction adiabatic multiplied by We\n ad3[k]:=add(Q[k,i]*la^ (i+k),i=0..P-k);\n s1:=s1+ad1[k]*d.k; s2:=s2+ad2[k]*q.k; s3:=s3+ad3[k ]*K.k;\nod;\n# sum of contributions adiabatic and nonadiabatic in the \+ energy\nDE:=simplify(We*(s1+s2)+s3); xtrct(P) end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "trunca:=proc(n) if n<0 then 0 else 1 fi: \+ end: # utility to truncate terms" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 564 "xtrct:=proc(P) # extraction of Z[k,l] \+ from energy\nlocal k,l,p,s,s1,z; global lista,Z,DE,vh,Be,We,JJ,g; \nfo r k from 0 to P/2 do s:=subs(vh=0,DE); DE:=expand((DE-s)/vh);\n for \+ l from 0 to P/2 do s1:=subs(JJ=0,s); s:=expand((s-s1)/JJ);\n for p f rom 0 to P do z[k,l,p]:=\n sort(collect(expand(subs(We=2*Be/g,coeff( s1,g,p)*g^p)),lista))\n od od od;\n for k from 0 to P/2 do\n f or l from 0 to P/2 do if k=0 and l=0 then Z[0,0]:=z[0,0,0]\n eli f z[k,l,k+2*l-1]<>0 and z[k,l,k+2*l]<>0\n then Z[k,l]:=z[k,l,k+2 *l-1]+z[k,l,k+2*l] fi od od end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 647 "convertacstu:=proc() # to convert from a_j to c_j, K_j to u_j etc.\nlocal j,l1,l2,l3,l4; global a,K,q,d,c,s,t,u,q0, t0,b0,s0;\nfor l1 from 1 to 12 do a.l1:=(-1)^l1*(l1+1)/2^l1+add((-1)^j *(l1+1)!\n *c.(l1-j)/(2^j*j!*(l1+1-j)!), j=0..l1 -1) od:\nfor l2 from 1 to 8 do K.l2:=add((-1)^(l2-j)*2^(j-l2)*(l2-1)! \n *u.j/((l2-j)!*(j-1)!),j=1..l2) od;\nfor l3 fr om 1 to 8 do q.l3:=add((-1)^(l3-j)*2^(j-l3)*(l3-1)!\n \+ *t.j/((l3-j)!*(j-1)!),j=1..l3) od; q0:=t0;\nfor l4 from 1 to 8 do d.l4:=add((-1)^(l4-j)*2^(j-l4)*(l4-1)!\n *s. j/((l4-j)!*(j-1)!),j=1..l4) od; b0:=s0; end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "tp3(6): Z[2,1]; # for instance" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 1 " " }{TEXT 261 31 "3 Vibrational matrix e lements" }}{PARA 0 "" 0 "" {TEXT -1 987 " Procedures in this suite VIBMATEL allow one to generate vibrational matrix elements of the coordinate x to various non-negative powers between vibrationa l states v and v'.\n To use this suite, specify the required order P of perturbation theory; then specify values of vibrational quantum \+ number for vibrational states v and v' as arguments of \"state\". For instance, to generate matrix elements of x to various powers between \+ vibrational states v=0 and v'=1 up to a_6, type \"P:=6: state(0): stat e(1):\". Particular matrix elements, in the array mel[v,l,v'], are fo rmed on typing \"me(v,l,v'):\". If matrix elements are req uired to\n be expressed in terms of coefficients c_j of function V(z) \+ for potential energy, after matrix elements are generated t hey can be so converted by invoking \"convertac():\". To ensure that \+ resulting expressions appear in attractive form, type \"series(sort(si mplify(mel[v',l,v])),g);\" to display a specific expression." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 85 "restart: P:=6: # Th e default value for greatest order of perturbation is 6. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 615 "state:=proc(n) # corrections t o energy e[n,p] and wavefunction c[j,n,p]\nlocal j,k,p,s,e: global c,P ,a:\nc[n,n,0]:=1:\nfor p from 0 to P do\nif p>0 then\n e[n,p]:=expand (add(a.s/2*add(Q(n,s+2,k)*c[k,n,p-s]*cut(k,n,p-s),\n k=max(n-s-2,0).. n+s+2),s=1..p)- add(e[n,s]*c[n,n,p-s],s=1..p-1))\nfi:\nfor j from 0 to n+3*p do\nif p>0 and j<>n then\n c[j,n,p]:=expand(1/(j-n)*(add(e[n,s ]*c[j,n,p-s]*cut(j,n,p-s),s=1..p)\n - add(a.s/2*add(c[k,n,p-s]*Q(j,s+ 2,k)*cut(k,n,p-s)\n ,k=max(j-s-2,0)..j+s+2),s=1..p))) fi: od:\nif p>0 then c[n,n,p]:=-expand(1/2*add(add(c[k,n,s]*c[k,n,p-s],\n k=0..n+3* min(s,p-s)),s=1..p-1)) fi: od: end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 223 "# matrix elements in a basis set of an har monic oscillator\nQ:=proc(i,n,j) option remember:\nif i<0 or j<0 then \+ 0\n elif n=0 then delta(i,j) else\n expand(sqrt(i/2)*Q(i-1,n-1,j)+sq rt((i+1)/2)*Q(i+1,n-1,j)) \nfi: end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 108 "# cutoff function to achieve c[j,n,p]=0 if |n-j|>3p \ncut:=proc(n,j,p) if abs(n-j)<=3*p then 1 else 0 fi: end:" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "delta:=proc(i,j) if \+ i=j then 1 else 0 fi end: # Kronecker delta" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 160 "me:=proc(j,k,n) global g,l,p,P,mel; # matrix elements \nl:=g^(1/2); mel[j,k,n]:=sum('mep(j,k,n,p)* l^p','p'=0..P):\nmel[n,k,j]:=mel[j,k,n]: end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 237 "# contributions to matrix elements \nm ep:=proc(m,k,n,p) local s,i,j: global c,l:\nl^k*sort(combine(expand(ad d(add(add(c[i,m,s]*c[j,n,p-s]*cut(i,m,s)\n *cut(j,n,p-s)*Q(i,k,j),j=0 ..n+3*(p-s)),i=0..m+3*s),s=0..p)),\n radical)): end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 194 "convertac:=proc() local j,k; globa l c; # to convert from a_j to c_j of V(z)\n for k from 1 to 12 do a .k:=(-1)^k*(k+1)/2^k+add((-1)^j*(k+1)!*\n c.(k-j)/(2^j*j!*(k+1-j) !), j=0..k-1) od end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "P :=6: state(0): state(1): me(0,1,1); # for instance" }{TEXT -1 0 "" }} }{EXCHG {PARA 0 "" 0 "" {TEXT 259 43 "4 Coefficients in the Herman-Wa llis factor" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 757 " Proc edures in this suite HERMWALL allow one to generate quantities used to form the Herman-Wallis factor F[0,v'](iota) for electric-dipolar tran sitions with J'-J= +1 or -1.\n To use this suite to form the Herma n-Wallis factor for the band from v=0 to v=v', invoke \"hw0(n):\" in w hich n is the numerical value of v'; then the results of the calculati on are alpha(0,v') and beta(0,v'), displayed on typing \"aln;\" and \" ben;\". For instance, for the band v=2 -- v=0, input \"hw0(2):\", the n \"al2;\" and \"be2;\". \+ \+ F[0,v'](iota) = 1 + 2 alv' i ota / <0|M(x)|v'> + (alv'^2 + 2 bev') iota^2 / <0|M(x)|v'>^2" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 789 "restart: hw0:=proc(v)\nloc al i,lista,dip: global g,l,JJ,m,P:\nP:=v+4: \+ # Greatest perturbation order\n# m=iota=1/2 [J'(J'+1)-J(J+1)] , J'=J+1,J-1 for electric-dipolar transitions\n# Matrix elements, ener gies etc. are reduced to terms in m^0, m^1, and m^2.\nlista:=[l,M.(0.. v)]:\nJJ:=m*(m-1): state(0):\nJJ:=m*(m+1): state(v):\n# first tw o terms of the matrix element <0J|M(x)|v'J'>\ndip:=l^(v+2)*add(M.i*mep (0,i,v,v+2-i),i=0..v+2)\n +l^(v+4)*add(M.i*mep(0,i,v,v+4-i),i=0..v+ 4):\ndip:=collect(dip,m):\ndip:=dip-coeff(dip,m,0):\n# Form alpha(0,v) and beta(0,v).\nal.v:=sort(collect(combine(coeff(dip,m,1)/l^v,radical ),lista)): al.0:=0:\nbe.v:=sort(collect(combine(coeff(dip,m,2)/l^v,rad ical),lista)):\nal.v:=subs(l=g^(1/2),al.v):\nbe.v:=subs(l=g^(1/2),be.v ): end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "delta:=proc(i, j) if i=j then 1 else 0 fi end: #Kronecker delta" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 233 "# matrix elements in a basis s et for an harmonic oscillator \nQ:=proc(i,n,j) option remember:\nif i< 0 or j<0 or n<0 then 0\n elif n=0 then delta(i,j) else\n expand(sqrt (i/2)*Q(i-1,n-1,j)+sqrt((i+1)/2)*Q(i+1,n-1,j)) fi: end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 112 "# truncation function to achieve c [j,n,p]=0 if |n-j|>3p\ncut:=proc(n,j,p) if abs(n-j)<=3*p then 1 else 0 fi: end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 99 "# truncation t o shift the the rotational part two orders\ncut2:=proc(s) if s>2 then \+ 1 else 0 fi end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1207 "# cor rections to energy e[n,p] and wavefunction c[j,n,p]\nstate:=proc(n) \n local j,k,p,s,e: global JJ,P,m,c:\nc[n,n,0]:=1:\nfor p from 0 to P do \n if p>0 then\n e[n,p]:=expand(add(a.s/2*add(Q(n,s+2,k)*c[k,n,p-s ]*cut(k,n,p-s),\n k=max(n-s-2,0)..n+s+2) + (-1)^s*(s-1)/2*JJ*cut2(s )*add(Q(n,s-2,k)\n *c[k,n,p-s]*cut(k,n,p-s),k=max(n-s,0)..n+s),s=1. .p)\n -add(e[n,s]*c[n,n,p-s],s=1..p-1)) fi:\n e[n,p]:=collect(e[ n,p],m):\n e[n,p]:=coeff(e[n,p],m,0)+m*coeff(e[n,p],m,1)+m^2*coeff(e[ n,p],m,2):\n for j from 0 to n+3*p do\n if p>0 and j<>n then\n c[ j,n,p]:=expand(1/(j-n)*(add(e[n,s]*c[j,n,p-s]\n *cut(j,n,p-s),s=1 ..p) - add(a.s/2*add(c[k,n,p-s]\n *Q(j,s+2,k)*cut(k,n,p-s),k=max( j-s-2,0)..j+s+2)\n +(-1)^s*(s-1)/2*JJ*cut2(s)*add(c[k,n,p-s]*Q(j, s-2,k)\n *cut(k,n,p-s),k=max(j-s,0)..j+s),s=1..p))) fi:\n c[ j,n,p]:=collect(c[j,n,p],m):\n c[j,n,p]:=coeff(c[j,n,p],m,0)+m*coeff( c[j,n,p],m,1)\n +m^2*coeff(c[j,n,p],m,2): od:\n \+ if p>0 then \n c[n,n,p]:=-expand(1/2*add(add(c[k,n,s]*c[k,n,p-s], \n k=0..n+3*min(s,p-s)),s=1..p-1)):\n c[n,n,p]:=collect(c[n,n,p] ,m):\n c[n,n,p]:=coeff(c[n,n,p],m,0)+m*coeff(c[n,n,p],m,1)\n +m^ 2*coeff(c[n,n,p],m,2):\n fi: od: end:" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 288 "# corrections to matrix elements \nmep:=p roc(mm,k,n,p) \nlocal s,i,j,s1: global c,m:\ns1:=expand(add(add(add(c[ i,mm,s]*c[j,n,p-s]*cut(i,mm,s)\n *cut(j,n,p-s)*Q(i,k,j),j=0..n+3*(p-s )),i=0..mm+3*s),s=0..p)):\ns1:=collect(s1,m):\ns1:=coeff(s1,m,0)+m*coe ff(s1,m,1)+m^2*coeff(s1,m,2):\nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 199 "convertac:=proc() local j,k; global c; # to c onvert from a_j to c_j of V(z)\n for k from 1 to 12 do a.k:=(-1)^k*(k +1)/2^k+add((-1)^j*(k+1)!*\n c.(k-j)/(2^j*j!*(k+1-j)!), j=0..k-1) od end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "hw0(2): al2; b e2; # for instance" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "Here endet h this worksheet." }}}}{MARK "34 0 0" 32 }{VIEWOPTS 1 1 0 1 1 1803 }