let u = [];

let f = [];

let F = [];

let r = [];

let P = [];

let p = [];

let x = [];

let y = [];

let R = [];

let Ri=[];

let charge= [];

let normu = [];

let N=200; 

let h=4/N;

let dt=0.05;

let M=5;

let M1=5;

let c=[];

let kinenergy = 0;

let potenergy = 0;

let totenergy = 0;

let D=0;

let t=0;

​

​

​

function setup() {

  createCanvas(400, 400);

  for (var m=1;m<M1;m++){

  w[m] = [];

  u[m] = [];

  P[m] = [];

  f[m] = [];

  c[m] = [];

  p[m] = [];

  r[m] = 0;

  x[m] = 0;

  y[m] = 0;

  R[m] =25;

  Ri[m] =0;

  D=20;

  charge[m]=1;

  normu[m] = 0;

  }

  for (var m=1;m<M1;m++){

  for (var i=0;i<N;i++){

  w[m][i]=[];

  u[m][i]=[];

  P[m][i]=[];

  p[m][i]=[];

  c[m][i]=[];

  f[m][i]=[];

  F[i]=[];

  for (var j=0;j<N;j++){

  w[m][i][j]=0;

  u[m][i][j]=0;

  P[m][i][j]=0;

  p[m][i][j]=0;

  c[m][i][j]=0;

  f[m][i][j]=0;

  F[i][j]=0;

  x[1]=100-2*D;

  y[1]=100-D;

  x[2]=100+2*D;

  y[2]=100-D;

  x[3]=100;

  y[3]=100+D;

  x[4]=100;

  y[4]=100+D;

  F[i][j]=-1/(h*sqrt(pow(i-x[1]+0.5,2)+pow(j-y[1],2))+4*h)-1/(h*sqrt(pow(i-x[2],2)+pow(j-y[2],2))+4*h)

    -1/(h*sqrt(pow(i-x[3],2)+pow(j-y[3],2))+4*h)

    -1/(h*sqrt(pow(i-x[4],2)+pow(j-y[4],2))+4*h)

  }

  }

  }

  Ri[3]=10;

  Ri[4]=10;

  for (var i=0;i<N;i++){

  for (var j=0;j<N;j++){

  for (var m=1;m<M;m++){

  r[m]=sqrt((i-x[m])*(i-x[m])+(j-y[m])*(j-y[m]));

  }

  if (r[1]<R[1]){

  w[1][i][j]=1;

  u[1][i][j]=1-r[1]/R[1];

  }

  if (r[2]<R[2]){

  w[2][i][j]=1;

  u[2][i][j]=1-r[2]/R[2];

  }

  if (r[3]<R[3] && j<100+D && r[3]>Ri[3]){

  w[3][i][j]=1;

  u[3][i][j]=1-r[3]/R[3];

  }

  if (r[4]<R[4] && j>100+D && r[4]>Ri[4]){

  w[4][i][j]=1;

  u[4][i][j]=1-r[4]/R[4];

  }

  }

  }

}

​

function draw() {

  background(255,255,255,255);

  if (t<1000){

  t=t+1;

​

  for (var i=1;i<N-1;i++){

  for (var j=1;j<N-1;j++){

  for (var m=1;m<M;m++){

  c[m][i][j]=0;

  for (var n=1;n<M;n++){

  if (n!=m){

   c[m][i][j]= c[m][i][j]-u[n][i][j];

  }

  }

  c[m][i][j]=0.5*(u[m][i][j]+c[m][i][j]);

  p[m][i][j]=0;

  for (var n=1;n<M;n++){

  if (n!=m){

  p[m][i][j]=p[m][i][j]+pow(u[n][i][j],2);

  }

  }

  }

  }

  }

  for (var q=1;q<10;q++){

  for (var m=1;m<M;m++){

  for (var i=1;i<N-1;i++){

  for (var j=1;j<N-1;j++){

  r[m]=sqrt((i-x[m])*(i-x[m])+(j-y[m])*(j-y[m]));

    if (r[m]>Ri[m]){

  w[m][i][j]=w[m][i][j]+0.5*dt*abs(c[m][i][j])*(w[m][i+1][j]+w[m][i-1][j]-4*w[m][i][j]+w[m][i][j+1]+w[m][i][j-1])+0.5*dt*c[m][i][j]*sqrt(pow((w[m][i+1][j]-w[m][i-1][j]),2)+pow((w[m][i][j+1]-w[m][i][j-1]),2)); 

 }

​

  u[m][i][j]=u[m][i][j]+0.5*dt*((u[m][i+1][j]-u[m][i][j])*(w[m][i+1][j]+w[m][i][j])/2-(u[m][i][j]-u[m][i-1][j])*(w[m][i][j]+w[m][i-1][j])/2+(u[m][i][j+1]-u[m][i][j])*(w[m][i][j]+w[m][i][j+1])/2-(u[m][i][j]-u[m][i][j-1])*(w[m][i][j]+w[m][i][j-1])/2)-dt*h*h*F[i][j]*w[m][i][j]*u[m][i][j]+dt*h*h*P[m][i][j]*w[m][i][j]*u[m][i][j];

P[m][i][j]=P[m][i][j]+0.05*((P[m][i+1][j]-4*P[m][i][j]+P[m][i-1][j]+P[m][i][j+1]+P[m][i][j-1])+0.05*h*h*200*p[m][i][j]);

  } 

  }

  }

  }

  }

for (var m=1;m<M;m++){

normu[m]=0;

for (var i=1;i<N-1;i++){

for (var j=1;j<N-1;j++){

if (w[m][i][j]>0.5){    

normu[m] = normu[m] + pow(u[m][i][j],2)*h*h;

}

}

}

for (var i=1;i<N-1;i++){

for (var j=1;j<N-1;j++){

u[m][i][j]= sqrt(charge[m])*u[m][i][j]/sqrt(normu[m]);

}

}

}

potenergy = 0;

kinenergy = 0;

for (var m=1;m<M;m++){

for (var i=1;i<N-1;i++){

for (var j=1;j<N-1;j++){

if (w[m][i][j]>0.5){    

potenergy = potenergy +(P[m][i][j]+F[i][j])*u[m][i][j]*u[m][i][j]*h*h

kinenergy = kinenergy + 0.5*(pow(u[m][i][j]-u[m][i+1][j],2)+pow(u[m][i][j]-u[m][i][j+1],2));

}

}

}

}

text(t,50,300);

text("time step",5,300);

text(potenergy,100,350);

text(kinenergy,100,370);

text(potenergy+kinenergy,100,390);

text("potential energy",10,350);

text("kinetic energy",10,370);

text("total energy",10,390);

for (var m=1;m<M;m++){

//ellipse(x[m]*2,y[m]*2,10);

}

​

noStroke();

for (var i=0;i<N;i++){

for (var j=0;j<N;j++){

for (var m=1;m<M;m++){

if (w[m][i][j]>0.5){ 

if (m==1){

fill(0,255,0,100*pow(u[1][i][j],2));

}

if (m==2){

fill(0,255,0,100*pow(u[2][i][j],2));

}

if (m==3){

fill(255,0,0,100*pow(u[3][i][j],2));

}

if (m==4){

fill(0,0,255,100*pow(u[4][i][j],2));

}

square(2*i,2*j,4);

}

}

}

}

​

for (var i=0;i<N;i++){

fill(0);

for (var m=1;m<M;m++){

fill(255,255,0,255);

circle(2*i,190-50*u[m][i][95],3);

fill(255,0,255,255);

circle(230-20*u[m][115][i],2*i,3);

fill(0,0,255,255);

circle(2*i,200-2*D-20*F[i][100+D],3);

fill(0,255,255,255);

circle(2*i,200-2*D+20*P[1][i][100-D],3);

fill(0);

circle(2*x[m],2*y[m],10);

​

}

}

​

  fill(0);

  textSize(10);

  text("RealQM 2d Molecule Model",10,20);

  text("Ansatz U(x)=sum_iU_i(x), non-overlapping supports",10,40);

  text("Minimize sum_i 0.5*integral(DU_i*DU_i-P*U_i) dx", 10,60);

  text("subject to integral U_i*U_i dx = charge_i ",10,80);

  text("Bernoulli Free Boundary: dU_i/dn =0,  U(x) continuous ",10,100);

  text("Blue: kernel potential, Lightblue: electron potential", 10, 120);

  text("Yellow: horisontal cross cut U(x), Black: support_i",10, 140);

  text("Magenta: vertical cross cut U(x)",10, 160);

​

​

}