#canned exercise for anharmonic 1d oscillator
 1    ;

function phi_i=harmonic_psi(x0,k,mass,i)
  #build parameters for harmonic solutions phi_i={x0,g,N,w,H}
  #phi_i(x)=N*e^{-(gx)^2/2}*H_i(gx)
  #N=sqrt(g/sqrt(pi)/2^i*i!)
  #w=sqrt(k/m)
  #g=(k*m)^1/4
  phi_i.i=i;
  phi_i.x0=x0;
  phi_i.w=sqrt(k/mass);
  phi_i.g=sqrt(sqrt(k*mass));
  phi_i.N=sqrt(phi_i.g/(pow2(i-1)*factorial(i-1)*sqrt(pi)));
  phi_i.H=hermite(i)*phi_i.N;
end  

function psi_x=psi_val(psi,x)
  xs=x*psi.g;
  psi_x=exp(-xs.*xs*0.5).*polyval(psi.H,xs);
end

function res=expectation(psi_a,psi_b,op)  
  xext=2*sqrt((psi_a.i+1)*(psi_b.i+1)/psi_a.g/psi_b.g); #tries to guess the real extent of the wf
  res=quadgk(@(x) op(x).*psi_val(psi_a,x).*psi_val(psi_b,x),psi_a.x0-xext,psi_a.x0+xext);
end

function y=op_id(x)
  y=1;
end

function y=op_x(x)
  y=x;
end

function y=op_x2(x)
  y=x.*x;
end

function [harmonic,anharmonic,v_harmonic,v_morse,e_harmonic,e_morse]=morse_solutions(e0,k,mass,n)
#returns n+1 analytic harmonic and anharmonic solutions in the morse potential of parameters e0,k
  omega=sqrt(k/mass);
  for i=1:n+1
    harmonic(i)=harmonic_psi(0,k,mass,i);
    e_harmonic(i)=omega*(i-0.5);
  endfor
  hm=zeros(n+1,n+1);
  v_harmonic=@(x) 0.5*k*x.*x;
  v_morse=morse_erk(e0,0,k,e0);
  for i=1:n+1
    for j=i:n+1
      thm=expectation(harmonic(i),harmonic(j),v_morse);
      thm=thm-expectation(harmonic(i),harmonic(j),v_harmonic);
      if i==j
        thm=thm+omega*(i-0.5);
      endif
      hm(i,j)=hm(j,i)=thm;
    endfor
  endfor
  # hm
  [q_morse,e_morse]=eigsort(hm);
  anharmonic=harmonic;
  for i=1:n+1
    H=[0];
    for j=1:n+1
      H=poly_add(H,harmonic(j).H*q_morse(j,i));
    endfor
    anharmonic(i).H=H*sign(H(1)*harmonic(i).H(1));
  endfor
endfunction

function plot_solutions(harmonic,anharmonic,v_harmonic,v_morse,e_harmonic,e_morse)
  ind_m=length(harmonic);
  psi_m=harmonic(ind_m);
  maxen=max(max(e_harmonic),max(e_morse))*4;
  xext=(psi_m.i+1)/psi_m.g/2; #tries to guess the real extent of the wf
  xx=linspace(-xext,xext,128);
  for i=1:ind_m
    y_harm=psi_val(harmonic(i),xx);
    y_morse=psi_val(anharmonic(i),xx);
    plot(xx,y_harm,xx,y_morse);
    set(gca,"xlim",[-xext,xext]);
    printf("la funzione d'onda di indice %d",i);
    pause();
  endfor
  rv_harm=v_harmonic(xx);
  rv_harm(rv_harm>maxen)=nan;
  rv_morse=v_morse(xx);
  rv_morse(rv_morse>maxen)=nan;
  plot(xx,y_harm);
  xxha=[-xext*0.5,0];
  xxmo=[0,xext*0.5];
  for i=1:ind_m
    l_harm(i,1:2)=e_harmonic(i);
    l_morse(i,1:2)=e_morse(i);
  endfor
  plot(xx,rv_harm,xx,rv_morse,xxha,l_harm,xxmo,l_morse);
  set(gca,"xlim",[-xext,xext]);
  input("il potenziale con i livelli energetici","s");
endfunction
