public class Complejo extends Object{
double real;
double imag;
static final Complejo i=new Complejo(0.0, 1.0);
public Complejo() {
real=0.0;
imag=0.0;
}
public Complejo(double real, double imag){
this.real=real;
this.imag=imag;
}
public Complejo opuesto(){
return new Complejo(-real, -imag);
}
public Complejo conjugado(){
return new Complejo(real, -imag);
}
public double modulo(){
return (real*real+imag*imag);
}
public double argumento(){
if(real==0.0 && imag==0.0){
return 0.0;
}else{
return Math.atan(imag/real);
}
}
public static Complejo suma(Complejo c1, Complejo c2){
double x=c1.real+c2.real;
double y=c1.imag+c2.imag;
return new Complejo(x, y);
}
public static Complejo diferencia(Complejo c1, Complejo c2){
double x=c1.real-c2.real;
double y=c1.imag-c2.imag;
return new Complejo(x, y);
}
public static Complejo producto(Complejo c1, Complejo c2){
double x=c1.real*c2.real-c1.imag*c2.imag;
double y=c1.real*c2.imag+c1.imag*c2.real;
return new Complejo(x, y);
}
public static Complejo cociente(Complejo c1, Complejo c2){
double aux=c2.real*c2.real+c2.imag*c2.imag;
double x=(c1.real*c2.real+c1.imag*c2.imag)/aux;
double y=(c1.imag*c2.real-c1.real*c2.imag)/aux;
return new Complejo(x, y);
}
public static Complejo cociente(Complejo c, double d){
double x=c.real/d;
double y=c.imag/d;
return new Complejo(x, y);
}
public static Complejo producto(Complejo c, double d){
double x=c.real*d;
double y=c.imag*d;
return new Complejo(x, y);
}
public static Complejo exponencial(Complejo c){
double x=Math.cos(c.imag)*Math.exp(c.real);
double y=Math.sin(c.imag)*Math.exp(c.real);
return new Complejo(x, y);
}
public static Complejo csqrt(double d){
if(d>=0) return new Complejo(Math.sqrt(d), 0);
return new Complejo(0, Math.sqrt(-d));
}
public String imprime(){
return new String(real+" "+imag+"*i\n");
}
}
public class CAux4x4 {
Complejo s[][]=new Complejo[2][4];
public CAux4x4() {
for(int i=0; i<2; i++){
for(int j=0; j<4; j++){
s[i][j]=new Complejo();
}
}
}
void auxiliar(Complejo q0, Complejo q1, double x){
s[0][0]=Complejo.exponencial
(Complejo.producto(Complejo.producto(Complejo.i, q0), x));
s[0][1]=Complejo.exponencial
(Complejo.producto(Complejo.producto(Complejo.i, q0), -x));
s[0][2]=Complejo.exponencial
(Complejo.producto(Complejo.producto(Complejo.i, q1), x)).opuesto();
s[0][3]=Complejo.exponencial
(Complejo.producto(Complejo.producto(Complejo.i, q1), -x)).opuesto();
s[1][0]=Complejo.producto(q0, s[0][0]);
s[1][1]=Complejo.producto(q0, s[0][1]).opuesto();
s[1][2]=Complejo.producto(q1, s[0][2]);
s[1][3]=Complejo.producto(q1, s[0][3]).opuesto();
}
}
public class CMatriz2x2 {
Complejo a11;
Complejo a12;
Complejo a21;
Complejo a22;
public CMatriz2x2() {
a11=new Complejo();
a12=new Complejo();
a21=new Complejo();
a22=new Complejo();
}
public CMatriz2x2(Complejo a11, Complejo a12, Complejo a21, Complejo a22){
this.a11=a11;
this.a12=a12;
this.a21=a21;
this.a22=a22;
}
public static CMatriz2x2 producto(CMatriz2x2 m1, CMatriz2x2 m2){
CMatriz2x2 temp=new CMatriz2x2();
temp.a11=Complejo.suma(Complejo.producto(m1.a11, m2.a11),
Complejo.producto(m1.a12,m2.a21));
temp.a12=Complejo.suma(Complejo.producto(m1.a11, m2.a12),
Complejo.producto(m1.a12,m2.a22));
temp.a21=Complejo.suma(Complejo.producto(m1.a21, m2.a11),
Complejo.producto(m1.a22,m2.a21));
temp.a22=Complejo.suma(Complejo.producto(m1.a21, m2.a12),
Complejo.producto(m1.a22,m2.a22));
return temp;
}
public static CVector2 producto(CMatriz2x2 m, CVector2 v){
CVector2 temp=new CVector2();
temp.a=Complejo.suma(Complejo.producto(m.a11, v.a),
Complejo.producto(m.a12,v.b));
temp.b=Complejo.suma(Complejo.producto(m.a21, v.a),
Complejo.producto(m.a22,v.b));
return temp;
}
public void matriz(Complejo k, double x){
a11=Complejo.exponencial
(Complejo.producto(Complejo.producto(Complejo.i, k), x));
a12=Complejo.exponencial
(Complejo.producto(Complejo.producto(Complejo.i, k), -x));
a21=Complejo.producto(k, a11);
a22=Complejo.producto(k.opuesto(), a12);
}
public void inversa(Complejo k, double x){
a11=Complejo.cociente(Complejo.exponencial(Complejo.producto
(Complejo.producto(Complejo.i, k), -x)), 2);
a21=Complejo.cociente(Complejo.exponencial(Complejo.producto
(Complejo.producto(Complejo.i, k), x)), 2);
a12=Complejo.cociente(a11, k);
a22=Complejo.cociente(a21, k.opuesto());
}
}
public class CVector2 {
Complejo a;
Complejo b;
public CVector2(Complejo a, Complejo b) {
this.a=a;
this.b=b;
}
public CVector2(){
a=new Complejo();
b=new Complejo();
}
public static CVector2 cociente(CVector2 v, double d){
return new CVector2(Complejo.cociente(v.a, d),
Complejo.cociente(v.b, d));
}
}
public class Matriz {
int filas;
int columnas;
Complejo a[][];
public Matriz(int filas, int columnas) {
this.filas=filas;
this.columnas=columnas;
a=new Complejo[filas][columnas];
for(int i=0; i<filas; i++){
for(int j=0; j<columnas; j++){
a[i][j]=new Complejo();
}
}
}
Complejo determinante(){
Complejo det=new Complejo(1.0, 0);
for(int k=0; k<filas-1; k++){
for(int i=k+1; i<filas; i++){
for(int j=k+1; j<filas; j++){
a[i][j]=Complejo.diferencia(a[i][j], Complejo.producto(a[i][k],
Complejo.cociente(a[k][j], a[k][k])));
}
}
det=Complejo.producto(det, a[k][k]);
}
det=Complejo.producto(det, a[filas-1][columnas-1]);
return det;
}
}
public abstract class Raices {
static final double E1=1.0e-15; //cota de errores
static final double E2=0.0001;
static final int MAXITERACCION=100;
static final int MAXNIVELES=30;
static final double pasos[]={0.2, 0.2, 0.2,
0.1, 0.1, 0.06, 0.05, 0.04, 0.03, 0.01};
int N; //número de barreras
double a; //anchura
double b; //separación
static final double H=5.0; //profundidad del pozo;
//resultados
double Ceros[]=new double[MAXNIVELES]; // ceros de la función
int nCeros=0; // número de ceros de la función
double DE; //paso de exploración en la búsqueda de ceros
abstract double f(double E, int paridad);
public Raices(int N, double a, double b) {
this.a=a;
this.b=b;
this.N=N;
DE=pasos[N];
}
//procedimiento del punto medio
double punto_medio(double inf, double sup, int paridad, Boolean bError){
double med, inter, ini;
int j=0; //cuenta el número de iteraciones
ini=f(inf, paridad);
do{
med=(inf+sup)/2;
inter=f(med, paridad);
if(Math.abs(inter)<E1) return med;
if(Math.abs((sup-inf)/med)<E2) return med;
if(inter*ini<0)
sup=med;
else{
inf=med;
ini=inter;
}
j++;
}while(!(j==MAXITERACCION));
bError=new Boolean(true);
return med;
}
//calcula en un intervalo entre dos discontinuidades
void calcular(double E_inf, double E_sup, int paridad){
if((E_inf>E_sup)||(Math.abs(E_sup-E_inf)<=DE)) return;
double inf, sup;
double ini, fin;
double raiz;
Boolean Error=new Boolean(false);
inf=E_inf+DE;
ini=f(inf, paridad);
sup=inf+DE;
while(sup<E_sup-DE){
fin=f(sup, paridad);
if(ini*fin>0){
ini=fin;
inf=sup;
sup=inf+DE;
continue;
}
raiz=punto_medio(inf, sup, paridad, Error);
if (!Error.booleanValue()){
Ceros[nCeros]=raiz;
nCeros++;
}
Error=new Boolean(false);
ini=fin;
inf=sup;
sup=inf+DE;
}
}
}
public class Niveles {
int paridad;
double energia;
public Niveles() {
paridad=0;
energia=0.0;
}
}
public class Solido extends Raices {
double x[]; //puntos de discontinuidad de la función potencial
double V[]; //función potencial
Complejo q[];
CVector2 Coef[]; //coeficientes para la función de onda
Matriz M; //matriz para el cálculo del determinante
int nNivel; //número de niveles totales
Niveles nivel[];
public Solido(int N, double a, double b) {
super(N, a, b);
V=new double[N+2];
M=new Matriz(2*N+2, 2*N+2);
//para la función de onda
q=new Complejo[N+2];
for(int j=0; j<N+2; j++){
q[j]=new Complejo();
}
Coef=new CVector2[N+2];
for(int j=0; j<N+2; j++){
Coef[j]=new CVector2();
}
//niveles de energía
nNivel=0;
nivel=new Niveles[MAXNIVELES];
for(int j=0; j<MAXNIVELES; j++){
nivel[j]=new Niveles();
}
//sólido regular
x=new double[N+3];
x[0]=0;
//número de barreras impar, número de átomos par
if(N%2==1){
x[1]=b/2;
for(int j=2;j<N+2; j+=2){
x[j]=x[j-1]+a;
if(j==N+1) break;
x[j+1]=x[j]+b;
}
for(int j=0; j<N+1; j+=2){
V[j]=H;
V[j+1]=0;
}
}else{
//número de barreras par, número de átomos impar
x[1]=a/2;
for(int j=2;j<N+2; j+=2){
x[j]=x[j-1]+b;
x[j+1]=x[j]+a;
}
for(int j=0; j<N+1; j+=2){
V[j]=0;
if(j==N) break;
V[j+1]=H;
}
}
V[N+1]=H;
//el 2 se puede cambiar significa hasta que la función
// de onda tiende a cero al hacerse x grande
x[N+2]=x[N+1]+2;
}
void inicializar(double E, int paridad){
CAux4x4 M4=new CAux4x4();
for(int j=0; j<N+2; j++){
q[j]=Complejo.csqrt(E-V[j]);
}
Complejo t1=new Complejo();
Complejo t2=new Complejo();
t1=Complejo.exponencial(Complejo.producto
(Complejo.producto(Complejo.i, q[0]), x[1]));
t2=Complejo.exponencial(Complejo.producto
(Complejo.producto(Complejo.i, q[0]), -x[1]));
M.a[0][0]=Complejo.suma(t1, Complejo.producto(t2, paridad));
M.a[1][0]=Complejo.producto(q[0],
Complejo.diferencia(t1, Complejo.producto(t2, paridad)));
t1=Complejo.exponencial(Complejo.producto
(Complejo.producto(Complejo.i, q[1]), x[1]));
M.a[0][1]=t1.opuesto();
M.a[1][1]=Complejo.producto(q[1].opuesto(), t1);
//si el número de barreras es mayor que cero
if(N>0){
t1=Complejo.exponencial(Complejo.producto
(Complejo.producto(Complejo.i, q[1]), -x[1]));
M.a[0][2]=t1.opuesto();
M.a[1][2]=Complejo.producto(q[1], t1);
t1=Complejo.exponencial(Complejo.producto
(Complejo.producto(Complejo.i, q[N]), x[N+1]));
M.a[2*N][2*N-1]=t1;
M.a[2*N+1][2*N-1]=Complejo.producto(q[N], t1);
t1=Complejo.exponencial(Complejo.producto
(Complejo.producto(Complejo.i, q[N]), -x[N+1]));
M.a[2*N][2*N]=t1;
M.a[2*N+1][2*N]=Complejo.producto(q[N].opuesto(), t1);
t1=Complejo.exponencial(Complejo.producto
(Complejo.producto(Complejo.i, q[N+1]), x[N+1])).opuesto();
M.a[2*N][2*N+1]=t1;
M.a[2*N+1][2*N+1]=Complejo.producto(q[N+1], t1);
}
for(int i=2; i<2*N; i+=2){
M4.auxiliar(q[i/2], q[i/2+1], x[i/2+1]);
for(int j=i; j<i+2; j++){
for(int k=i; k<i+4; k++){
M.a[j][k-1]=M4.s[j-i][k-i];
}
}
}
}
//define la función abstracta
double f(double E, int paridad){
Complejo res=new Complejo();
inicializar(E, paridad);
res=M.determinante();
//la parte imaginaria es pequeña
return (res.real+res.imag);
}
//método de la burbuja
//ordena según el valor de la energía
void ordenar(){
double d;
int j;
for(int izq=1; izq<nNivel; ++izq){
for(int der=nNivel-1; der>=izq; --der){
if(nivel[der-1].energia>nivel[der].energia){
d=nivel[der-1].energia;
j=nivel[der-1].paridad;
nivel[der-1].energia=nivel[der].energia;
nivel[der-1].paridad=nivel[der].paridad;
nivel[der].energia=d;
nivel[der].paridad=j;
}
}
}
}
//halla todos los niveles de energía
void hallarNiveles(){
nNivel=0;
for(int par=-1; par<=1; par+=2){
nCeros=0;
calcular(0.0, H, par);
for(int i=0; i<nCeros; i++){
nivel[nNivel+i].energia=Ceros[i];
nivel[nNivel+i].paridad=par;
}
nNivel+=nCeros;
}
ordenar();
}
boolean normalizar(double E, int paridad){
CVector2 X=new CVector2(new Complejo(1.0,0), new Complejo(paridad, 0));
CVector2 Y=new CVector2(new Complejo(1.0,0), new Complejo());
Coef[0]=X;
Coef[N+1]=Y;
CMatriz2x2 R=new CMatriz2x2();
CMatriz2x2 I=new CMatriz2x2();
for(int j=0; j<N+2; j++){
q[j]=Complejo.csqrt(E-V[j]);
//no se puede normalizar la función de onda
if(q[j].modulo()<1.0e-20) return false;
}
for(int j=1; j<N+2; j++){
I.inversa(q[j], x[j]);
R.matriz(q[j-1], x[j]);
X=CMatriz2x2.producto(I, CMatriz2x2.producto(R, X));
Coef[j]=X;
}
Complejo norma=new Complejo();
Complejo temp1=new Complejo();
Complejo temp2=new Complejo();
for(int j=0; j<N+1; j++){
if(E<V[j]){
//a y b son números reales y q[j] imaginario
norma=Complejo.suma(norma, Complejo.producto
(Complejo.producto(Coef[j].a, Coef[j].b), 2*(x[j+1]-x[j])));
temp1=Complejo.cociente(Complejo.producto
(Coef[j].a, Coef[j].a.conjugado()), Complejo.producto
(Complejo.producto(Complejo.i, q[j]), 2));
temp2=Complejo.diferencia(Complejo.exponencial
(Complejo.producto(Complejo.producto(Complejo.i, q[j]), 2*x[j+1])),
Complejo.exponencial(Complejo.producto
(Complejo.producto(Complejo.i, q[j]), 2*x[j])));
norma=Complejo.suma(norma, Complejo.producto(temp1, temp2));
temp1=Complejo.cociente(Complejo.producto
(Coef[j].b, Coef[j].b.conjugado()), Complejo.producto
(Complejo.producto(Complejo.i, q[j]), 2));
temp2=Complejo.diferencia(Complejo.exponencial
(Complejo.producto(Complejo.producto(Complejo.i, q[j]), -2*x[j+1])),
Complejo.exponencial(Complejo.producto
(Complejo.producto(Complejo.i, q[j]), -2*x[j])));
norma=Complejo.diferencia(norma, Complejo.producto(temp1, temp2));
}else{
//a y b son complejos conjugados y q[j] es real
temp1=Complejo.suma(Complejo.producto(Coef[j].a,
Coef[j].a.conjugado()), Complejo.producto(Coef[j].b, Coef[j].b.conjugado()));
norma=Complejo.suma(norma, Complejo.producto(temp1, 2*(x[j+1]-x[j])));
temp1=Complejo.cociente(Complejo.producto(Coef[j].a,
Coef[j].b.conjugado()), Complejo.producto(Complejo.producto(Complejo.i, q[j]), 2));
temp2=Complejo.diferencia(Complejo.exponencial
(Complejo.producto(Complejo.producto(Complejo.i, q[j]), 2*x[j+1])),
Complejo.exponencial(Complejo.producto
(Complejo.producto(Complejo.i, q[j]), 2*x[j])));
norma=Complejo.suma(norma, Complejo.producto(temp1, temp2));
temp1=Complejo.cociente(Complejo.producto(Coef[j].b,
Coef[j].a.conjugado()), Complejo.producto(Complejo.producto(Complejo.i, q[j]), 2));
temp2=Complejo.diferencia(Complejo.exponencial
(Complejo.producto(Complejo.producto(Complejo.i, q[j]), -2*x[j+1])),
Complejo.exponencial(Complejo.producto
(Complejo.producto(Complejo.i, q[j]), -2*x[j])));
norma=Complejo.diferencia(norma, Complejo.producto(temp1, temp2));
}
}
temp1=Complejo.cociente(Complejo.producto(Coef[N+1].a,
Coef[N+1].a.conjugado()), Complejo.producto(Complejo.producto(Complejo.i, q[N+1]), -2));
temp2=Complejo.exponencial(Complejo.producto
(Complejo.producto(Complejo.i, q[N+1]), 2*x[N+1]));
norma=Complejo.suma(norma, Complejo.producto(temp1, temp2));
for(int j=0; j<N+2; j++){
Coef[j]=CVector2.cociente(Coef[j], Math.sqrt(norma.real));
}
return true;
}
}
|
