# include <stdlib.h>
# include <stdio.h>
# include <time.h>

//# include "eispack.h"
#include "eispack.c"
//Resuelve la ec de Schrodinger unidimensional. Particula confinada.
//iarmonico !=0 osc armónico, omega = frecuencia angular
//iarmonico=0 pozo cuadrado infinito de  potencial (a, b) en metros: paredes

//variables globales 
char iarmonico=0;
int nf=5; //funcion de onda para listar, 0 = fundamental
double h = 1.05459e-34; //h barra
double m = 1./6.02217e26; // masa en kg, numerador en uma 
double a=0.e-10,b=1.e-10, omega=1.0e14; //limites del pozo y w del armonico
double V(double x)
// Funcion potencial. Se toma el mínimo como cero
// V<0 => infinito
//retorna 2mV(x)/h^2
{    
      //output: 2*m*V/h^2;
   if(iarmonico) return m*omega/h*x*m*omega/h*x;  //Osc Armonico

   else{
      if(x>=a&&x<=b) return 0.;// Pozo cuadrado de altura infinita
      else return -1.; // Retorna negativo para significar infinito
/* Anadir la definicion para otro potencial.*/
      }
}       
main()
{/* Metodo de diferencias finitas.
    Sea la funcion de onda F(x)
    y la Ec. de Schrodinger 1d: [-h^2/m d^2/dx^2-V(x)]F(x)=EF(x)
    Se aproxima  d^2F/dx^2 = (Fi+1+Fi-1-2*Fi)/Dx^2
    Definiendo alpha(x) = (2*m*Dx^2/h^2)*V(x) y lambda = (2*m*Dx^2/h^2)*E
    Se cumple la ecuacion en diferencias finitas:
    (2+alpha(x_i))*Fi -Fi+1-Fi-1 =lambda*Fi
    Considerando las Fi como componentes de un vector, 
    es una ecuación de autovalores y autovectores.
    
    Se resuelve numéricamente usando el paquete libre eispack (de GNU) que
    obtiene los autovalores y autovectores (valores de la autofuncion para cada
    autovalor)
    Se supone que el potencial debe ser muy alto excepto para x1 < x < x2.
    unico rango de valores con F ditinta de cero. 
  */  
FILE * pfileout;
 int n=1001,err,i,j,ilist=51; double *d,*e,*z;
 double pi=3.141592653589793238;
 double Dx; // espaciado en metros entre puntos para calcular la funciond e onda
 double EeV,Ehw,EdivE0;
 double zteo;
 double x1,x2,xi; //limites de x para calcular
 double Norma;
 pfileout=fopen("Schrodinger.txt","w");
 if(iarmonico==1)
 {
 //para  Osc Arm, hacemos Fi=0 si  V > 100*h*omega
 x2=sqrt(h/m/omega)*sqrt(200.); 
 x1=-x2;
 Dx=(x2-x1)/n;
 }
 else
 {
 Dx=(b-a)/(n+1);
 x1=a+Dx;
 x2=b-Dx;
  omega=h/m/(b-a)/(b-a)*pi*pi/2.; //omega = E0/h teórico, para pozo cuadrado
  //realmente omega es un valor irrelevante
 }
 printf("\n Datos del calculo:");
 printf("\n Tipo de potencial: ");
 fprintf(pfileout,"\n Datos del calculo:");
 fprintf(pfileout,"\n Tipo de potencial: ");
if(iarmonico){printf(" OSCILADOR ARMONICO");fprintf(pfileout," OSCILADOR ARMONICO");}
 else {printf(" POZO CUADRADO");fprintf(pfileout," POZO CUADRADO");}
 printf("\n m =%10le = %5.3lf uma,h = %le (hbarra)",m,m*6.02217e26,h);
 fprintf(pfileout,"\n m =%10le = %5.3lf uma",m,m*6.02217e26);
 printf("\n omega: %le rad/s = E1(anal)/h para pozo cuadrado",omega);
 fprintf(pfileout,"\n omega: %le rad/s = E1(anal)/h para pozo cuadrado",omega);
 if(iarmonico){
 printf("\n omega: %10le rad/s",omega);
 printf("\n limites de x para calculo: %le %le metros",x1,x2);
 fprintf(pfileout,"\n omega: %10le rad/s",omega);
 fprintf(pfileout,"\n limites de x para calculo: %le %le metros",x1,x2);
 }
 else {
      printf("\n Paredes del pozo cuadrado: %le %le metros",a,b);
      fprintf(pfileout,"\n Paredes del pozo cuadrado: %le %le metros",a,b);    
      }
 printf("\n Dx: %le m, no. de puntos calculados de Fi: %i",Dx,n); 
 fprintf(pfileout,"\n Dx: %le m, no. de puntos calculados de Fi: %i",Dx,n);
 printf("\nteclear <RET>");
 getchar();//para que no se cierre la ventana
 /*
 d = entrada: diagonal mayor de la matriz 2+alpha[i]
       salida: vector de autovalores
 e = diagonales superior e inferior e[i]=-1, e[0] no se usa
 z= matriz Para convertir matriz simétrica en tridiagonal. 
  Para tridiagonal (es el caso) z no se usa y se da la matriz identidad listada por filas.
*/
// 
 d = ( double * ) malloc ( n* sizeof ( double ) );
 e = ( double * ) malloc ( n* sizeof ( double ) );
 z = ( double * ) malloc ( n*n* sizeof ( double ) );
 /* ejemplo diagonalizacion
 for(i=0;i<n;i++)
 {
 for(j=0;j<n;j++)z[i+n*j]=0.;
 d[i]=2.; 
 e[i]=-1.;
 z[i+n*i]=1.;
  }
  */
  printf("\n Lista algunos valores de la diagonal:2+ 2*m*Vi*Dx^2/h^2");
  printf("\nNum  2+alpha     x[i](m)     2mV(xi)/h**2      eV");
  printf("\n=================================================");
  fprintf(pfileout,"\n Lista algunos valores de la diagonal:2+ 2*m*Vi*Dx^2/h^2");
  fprintf(pfileout,"\nNum  2+alpha     x[i](m)     2mV(xi)/h**2      eV");
  fprintf(pfileout,"\n=================================================");
 for(i=0;i<n;i++)
 {
 xi =x1+i*Dx; 
 for(j=0;j<n;j++)z[i+n*j]=0.;
 d[i]=2.+V(xi)*Dx*Dx; 
 e[i]=-1.;
 /* z[] matriz identidad n x n. Listada por filas
  En este programa sólo  lista los autovalores y un autovector
 */
 z[i+n*i]=1.; 
 EeV=V(xi)*h*h/2/m/1.60217e-19;
 if(!(i%((n-1)/10))) //Lista solo unos cuantos valores de las Vi
 {
 printf("\n %i %10lf %10le %le %lf",i,d[i],xi,V(xi),EeV);
 fprintf(pfileout,"\n %i %10lf %10le %le %lf",i,d[i],xi,V(xi),EeV);
 }
/* if(i%100==0) //para listar con pausa cada 100 lineas
   {
    printf(" teclear <RET>");
    getchar();
   }
*/
  }
 printf("\n\n Diagonalizando...");
 
 err=tql2(n,d,e,z);
 printf("diagonalizacion finalizada");
 printf("\n\n Lista %i primeros autovalores y energias",ilist);
 printf("\nn  lambda_n    En(eV)     En/hw     E/E0");
 fprintf(pfileout,"\n Lista %i primeros autovalores y energias",ilist);
 fprintf(pfileout,"\nn  lambda_n    En(eV)     En/hw     E/E0");
 for(i=0;i<ilist;i++) 
 {
 EeV=d[i]/2*h/m/Dx*h/Dx/1.60217e-19;
 Ehw=EeV*1.60217e-19/h/omega;
 EdivE0=d[i]/d[0];
 printf("\n%i %10lf %10lf %10lf %10lf",i,d[i],EeV,Ehw,EdivE0);
 fprintf(pfileout,"\n%i %10lf %10lf %10lf %10lf",i,d[i],EeV,Ehw,EdivE0); }

 printf("\n\n\n Teclear <RET> para listar funcion de onda No %i",nf);
 getchar();
 printf("\n\n no punto Fi(x_I)       Fi_analitica");
 fprintf(pfileout,"\n\n no punto Fi(x_I)       Fi_analitica"); 
 //Normaliza funcion no nf: integral trapecio
 Norma=0;
 for(i=0;i<n;i++) Norma=Norma+z[i+nf*n]*z[i+nf*n];
 Norma=sqrt(Dx*(Norma-(z[0+nf*n]*z[0+nf*n]+z[n-1+nf*n]*z[n-1+nf*n])/2.));
 
 for(i=0;i<n;i+=(n-1)/100)
 {
 // Lista funcion de onda analitica. Armonico solo fundamental
 if(iarmonico)zteo = 
  sqrt(sqrt(m*omega/h/pi))*exp(-m*omega/h*(x1+i*Dx)*(x1+i*Dx)/2);
  // Pozo cuadrado funcion de onda analitica para nivel nf
 else zteo=sqrt(2./(b-a))*sin((nf+1)*pi*((i+1)*Dx)/(b-a)); 
 printf("\n %5i\t%10.2lf\t%10.2lf",i,z[i+nf*n]/Norma,zteo);
 fprintf(pfileout,"\n %5i\t%10.2lf\t%10.2lf",i,z[i+nf*n]/Norma,zteo);
 } 
  
 free ( d );
 free ( e );
 free ( z );
 printf("\n Teclear <RET> para terminar");
 getchar();
 fclose(pfileout); 
}