createCanvas(400, 400);

}

​

function draw() {

  background(220);

}

​

let u = [];

let K = [];

let Z;

let E1;

let E1red;

let E2;

let E3;

let E2red;

let E3red;

let N=400;

let M0=40;

let M1=100;

let M2=200;

let dt=1;

let h;

let normu1=0;

let normu2=0;

let normu3=0;

let energy1=0;

let energy2=0;

let energy3=0;

let slider1;

let slider2;

let slider3;

let slider0;

let start;

let stop;

let P = [];

let fix=0;

let D=2;

let S=0;

let diff=0.001;

​

function setup() {

  createCanvas(400, 400);

  slider1=createSlider(2,2,2,1);

  slider2=createSlider(0,8,8,1);

  slider3=createSlider(0,18,18,1);

  slider0=createSlider(2,28,28,1);

  start=createButton('Start');

  stop=createButton('Stop');

  start.position(350,365);

  start.mousePressed(startsimulation);

  stop.mousePressed(stopsimulation);

  stop.position(350,385);

  h=D/N;

  Z=slider0.value();

  E1=slider1.value();

  E2=slider2.value();

  E3=slider3.value();

  E1red=(E1-1)/E1;

  E2red=(E2-1)/E2;

  E3red=(E3-1)/E3;

  dt=0.5*pow(h,2);

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

  u[i]=exp(-Z*(i*h));

  P[i]=0;

  }

  u[N]=0;

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

  K[i]=Z/(i*h+fix*h);

  }

  for (var i=M0; i<M1+1;i++){

  if (i>1){

  for (var j=1;j<i;j++){       P[i]=P[i]+0.5*E1red*pow(u[j],2)*pow(j*h,2)*h/(i*h);

  }

  }

  for (var j=i;j<M1+1;j++){

  P[i]=P[i]+0.5*E1red*pow(u[j],2)*pow(j*h,2)*h/(j*h);

  }

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

  P[i]=P[i]+0.5*pow(u[j],2)*pow(j*h,2)*h/(j*h);

  }

  }

  for (var i=M1+1; i<M2+1;i++){

  P[i]=0;

  if (i>M1+1){ 

  for (var j=M1+1;j<i;j++){       P[i]=P[i]+0.5*E2red*pow(u[j],2)*pow(j*h,2)*h/(i*h);

  }

  }

  for (var j=i;j<M2+1;j++){

  P[i]=P[i]+0.5*E2red*pow(u[j],2)*pow(j*h,2)*h/(j*h);

  }

  for (var j=1;j<M1+1;j++){       P[i]=P[i]+0.5*pow(u[j],2)*pow(j*h,2)*h/(i*h);

  }

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

  P[i]=P[i]+0.5*pow(u[j],2)*pow(j*h,2)*h/(j*h);

  }

  }

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

  P[i]=0;

  if (i>M2+1){ 

  for (var j=M2+1;j<i;j++){       P[i]=P[i]+0.5*E2red*pow(u[j],2)*pow(j*h,2)*h/(i*h);

  }

  }

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

  P[i]=P[i]+0.5*E2red*pow(u[j],2)*pow(j*h,2)*h/(j*h);

  }

  for (var j=1;j<M2+1;j++){       P[i]=P[i]+0.5*pow(u[j],2)*pow(j*h,2)*h/(i*h);

  }

  }

}

​

​

function draw() {

  background(220);

  if (S==0){

  fill(0,0,255,100);

  circle(200,200,300);

  fill(0,255,0,255);

  circle(200,200,150);

  fill(255,0,0,255);

  circle(200,200,50);

  fill(0);

  circle(200,200,5);

  Z=slider0.value();

  E1=slider1.value();

  E2=slider2.value();

  E3=slider3.value();

  noStroke();

  fill(0);

  text("E1",205,200);

  text("E2",240,200);

  text("E3", 290,200);

  text("RealQM Atom Simulator: Ground States + 1st,2nd Ions",50,20);

  text("Kernel Charge 10 < Z < 19.",100,40);

  text("1st Shell: 2, 2nd Shell: E2, 3rd Shell: E3 Electrons.",100,60);

  text("Choose Z, E2 and E3 below and press Start!",100,130);

  text("RealQM = System of nonoverlapping electron densities",100,80);

  text("meeting at Bernoulli free boundary.",160,100);

  text("Electronic shell struture in accordance with observation:",50,230);

  text("Lithium Z=3: 2+1 (7.48),1st 2 (7.28)",50,250);

  text("Beryllium Z=4: 2+2 (14.6), 1st 2+1 (14.2), 2nd 2 (13.5)",50,265);

  text("Boron Z=5: 2+2+1 (24.5),1st 2+2 (24.2)",50,280);

  text("Carbon Z=6: 2+2+2 (37.5),1st 2+2+1 (37.1), 2nd 2+2 (36.5)",50,295);

  text("Nitrogen Z=7: 2+3+2 (54.6), 1st 2+2+2 (54.1)",50,310);

  text("Oxygen Z=8: 2+3+3 (74.8), 1st 2+3+2 (74.2)",50,325);

  text("Fluorine Z=9: 2+4+3 (99.8), 1st 2+4+2 (99.2)",50,340);

  text("Neon Z=10: 2+4+4 (128.5), 1st 2+4+3 (127.7), 2nd 2+4+2 (126.2)",50,355);

  text("Z =",50,399);

  text(Z,70,399);

  text("E2 =",180,399);

  text("E3 =",280,399);

  text(E2,210,399);

  text(E3,310,399);

  if (E1==0){

  M0=400;

  }

  if (E2==0){

  M1=400;

  }

  if (E3==0){

  M2=400;

  }

  if (Z==3 && E2==1 && E3==0){

//  text("(Li obs = -7.48)",285,365);

  fix=0;

  D=5;

  }

  if (Z==3 && E2==0){

//  text("(Li+ obs = -7.28)",285,365);

  fix=0;

  D=7;

  }

  if (Z==4 && E2==2){

//  text("(Be obs = -14.5)",280,365);

  fix=0;

  D=7;

  }

  if (Z==4 && E2==1){

//  text("(Be+ obs = -14.2)",280,365);

  fix=0;

  D=7;

  }

  if (Z==4 && E2==0){

//  text("(Be2+ obs= -13.5)",270,365);

  fix=0;

  D=7;

  }

 if (Z==5 && E2==2 && E3==1){

//  text("(Bo obs= -24.53)",285,365);

  fix=0.1

  D=7;

  }

  if (Z==5 && E2==2 && E3==0){

//  text("(Bo+ obs= -24.53)",285,365);

  fix=0.7;

  D=7;

  }

  if (Z==6 && E2==2 && E3==2){

//  text("(C obs= -37.7)",285,365);

  fix=0;

  D=7;

​

  }

  if (Z==7 && E2==3 && E3==2){

//  text("(N obs= -54.6)",285,365);

  fix=0;

  D=7;

  }

  if (Z==7 && E2==4 && E3==0){

//  text("(N+ obs= -54.1)",285,365);

  fix=0.60;

  D=7;

  }

  if (Z==8 && E2==3 && E3==3){

//  text("(O obs = -74.8)",285,365);

  fix=0;

  D=7;

  }

  if (Z==8 && E2==3 && E3==2){

//  text("(O+ obs = -74.2)",280,365);

  fix=0;

  D=7;

  }

  if (Z==8 && E2==2 && E3==2){

//  text("(O+ obs = -73.0)",270,365);

  fix=0;

  D=7;

  }

  if (Z==9 && E2==4 && E3==3){

//  text("(Fluor obs = -99.8)",270,365);

  fix=0.07;

  D=7;

  }

  if (Z==9 && E2==4 && E3==3){

//  text("(Fluor+ obs = -99.1)",270,365);

  fix=0.07;

  D=7;

  }

  if (Z==10 && E2==4 && E3==4){

//  text("(Neon obs = -128.5)",270,365);

  fix=0.1;

  D=7;

  }

  if (Z==10 && E2==4 && E3==3){

// text("(Neon+ obs = -127.7)",270,365);

  fix=0.1;

  D=7;

  }

  if (Z==10 && E2==4 && E3==2){

// text("(Neon2+ obs = -126.2)",270,365);

  fix=0.1;

  D=7;

  }

  h=D/N;

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

  K[i]=Z/(i*h+fix*h);

  }

  }

​

if (S==1){ 

  if (u[M1]>u[M1+1]+diff){  

   M1=M1+1;

  }

  if (u[M1]<u[M1+1]-diff){

  M1=M1-1;

  }

  if (u[M2]>u[M2+1]+diff){  

   M2=M2+1;

  }

  if (u[M2]<u[M2+1]-diff){

  M2=M2-1;

  }

​

for (m=0;m<20;m++){

for (k=0;k<5;k++){

  u[0]=u[1]/(1-Z*h);

  energy1=0;

  u[M1]=u[M1-1];

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

u[i]=u[i]+0.5*dt*(u[i+1]-2*u[i]+u[i-1])/pow(h,2)+0.5*dt*(u[i+1]-u[i-1])/(h*i*h) + dt*K[i]*u[i] - dt*2*P[i]*u[i];

  energy1= energy1 + 0.5*pow((u[i+1]-u[i -1])/(2*h),2)*pow(i*h,2)*h - pow(u[i],2)*K[i]*pow((i*h),2)*h + P[i]*pow(u[i],2)*pow(i*h,2)*h;

   }

  energy2=0;

  u[M1+1]=u[M1+2];

  u[M2]=u[M2-1];

   for (var i=M1+2;i<M2;i++){

   u[i]=u[i]+0.5*dt*(u[i+1]-2*u[i]+u[i-1])/pow(h,2)+0.5*dt*(u[i+1]-u[i])/(h*i*h) + dt*K[i]*u[i] - dt*2*P[i]*u[i];

  energy2= energy2 + 0.5*pow((u[i+1]-u[i-1])/(2*h),2)*pow(i*h,2)*h -pow(u[i],2)*K[i]*pow((i*h),2)*h + P[i]*pow(u[i],2)*pow(i*h,2)*h;

   }

  energy3=0;

  u[M2+1]=u[M2+2];

   for (var i=M2+2;i<N;i++){

   u[i]=u[i]+0.5*dt*(u[i+1]-2*u[i]+u[i-1])/pow(h,2)+0.5*dt*(u[i+1]-  u[i])/(h*i*h) + dt*K[i]*u[i] - dt*2*P[i]*u[i];

  energy3= energy3 + 0.5*pow((u[i+1]-u[i-1])/(2*h),2)*pow(i*h,2)*h -pow(u[i],2)*K[i]*pow((i*h),2)*h + P[i]*pow(u[i],2)*pow(i*h,2)*h;

   }

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

  P[i]=0;

  if (i>1){

  for (var j=1;j<i;j++){           P[i]=P[i]+0.5*E1red*pow(u[j],2)*pow(j*h,2)*h/(i*h);

  }

  }

  for (var j=i;j<M1+1;j++){

  P[i]=P[i]+0.5*E1red*pow(u[j],2)*pow(j*h,2)*h/(j*h);

  }

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

  P[i]=P[i]+0.5*pow(u[j],2)*pow(j*h,2)*h/(j*h);

  }

  }

  for (var i=M1+1; i<M2+1;i++){

  P[i]=0;

  if (i>M1+1){ 

  for (var j=M1+1;j<i;j++){       P[i]=P[i]+0.5*E2red*pow(u[j],2)*pow(j*h,2)*h/(i*h);

  }

  }

  for (var j=i;j<M2+1;j++){

  P[i]=P[i]+0.5*E2red*pow(u[j],2)*pow(j*h,2)*h/(j*h);

  }

  for (var j=1;j<M1+1;j++){       P[i]=P[i]+0.5*pow(u[j],2)*pow(j*h,2)*h/(i*h);

  }

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

  P[i]=P[i]+0.5*pow(u[j],2)*pow(j*h,2)*h/(j*h);

  }

  }

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

  P[i]=0;

  if (i>M2+1){ 

  for (var j=M2+1;j<i;j++){       P[i]=P[i]+0.5*E2red*pow(u[j],2)*pow(j*h,2)*h/(i*h);

  }

  }

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

  P[i]=P[i]+0.5*E2red*pow(u[j],2)*pow(j*h,2)*h/(j*h);

  }

  for (var j=1;j<M2+1;j++){       P[i]=P[i]+0.5*pow(u[j],2)*pow(j*h,2)*h/(i*h);

  }

  }

}

  normu1=0;

  for (var i=0;i<M1+1;i++){ 

  normu1= normu1 + pow(u[i]*(i*h),2)*h;

  }

  for (var i=0;i<M1+1;i++){ 

  u[i]=sqrt(E1)*u[i]/sqrt(normu1);

  }

  normu2=0;

  for (var i=M1+1;i<M2+1;i++){ 

  normu2= normu2 + pow(u[i]*(i*h),2)*h;

  }

  for (var i=M1+1;i<M2+1;i++){ 

   u[i]=sqrt(E2)*u[i]/sqrt(normu2);

  }

  normu3=0;

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

  normu3= normu3 + pow(u[i]*(i*h),2)*h;

  }

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

   u[i]=sqrt(E3)*u[i]/sqrt(normu3);

  }

}

  noStroke();

  if (Z==3 && E2==1){

  text("(Li obs= 7.48)",285,365);

  }

  if (Z==3 && E2==0){

  text("(Li+ obs = 7.28)",285,365);

  }

  if (Z==4 && E2==2){

  text("(Be obs = 14.5)",280,365);

  }

  if (Z==4 && E2==1){

  text("(Be+ obs = 14.2)",280,365);

  }

  if (Z==4 && E2==0){

  text("(Be2+ obs = 13.5)",280,365);

  }

  if (Z==5 && E2==2 && E3==1){

  text("(Bo obs = 24.5)",285,365);

  }

  if (Z==5 && E2==0){

  text("(Bo+ obs = 24.2)",285,365);

  }

  if (Z==6 && E2==2 && E3==2){

  text("(C obs = 37.7)",285,365);

  }

  if (Z==7 && E2==3 && E3==2){

  text("(N observed= 54.6)",285,365);

  }

  if (Z==7 && E2==2 && E3==2){

  text("(N+ observed= 54.1)",285,365);

  }

  if (Z==8 && E2==3 && E3==3){

  text("(O observed = 74.8)",285,365);

  }

  if (Z==8 && E2==3 && E3==2){

  text("(O+ observed = 74.3)",280,365);

  }

  if (Z==8 && E2==2 && E3==2){

  text("(O2+ observed = 73.0)",280,365);

  }

  if (Z==9 && E2==4 && E3==3){

  text("(Fluor = 99.8)",285,365);

  }

  if (Z==9 && E2==4 && E3==2){

  text("(Fluor+ = 99.2)",285,365);

  }

  if (Z==10 && E2==4 && E3==4){

  text("(Neon observed = 128.5)",275,365);

  fix=0;

  D=7;

  }

  if (Z==10 && E2==4 && E3==3){

  text("(Neon+ = 127.7)",275,365);

  fix=0

  D=7;

  }

  if (Z==10 && E2==4 && E3==2){

  text("(Neon2+ = 126.2)",275,365);

  fix=0

  D=7;

  }

  noStroke();

  fill(0);

  text("3(sub)Shell Atoms: Ground State+Ions SpherSym.",50,10);

  text("Kernel K(x)= Z/|x|, 2<Z<11, u(r) wave function,",50,30);

  text("u(r)^2 density of electrons in shell1, shell2,3 subshells, r=|x|.",20,50);

  text("Normalisation: 0<E2<9, 0<=E3<5.",50,70);

   text("integral_1 u(r)^2*r^2*dr = 2, integral_2,3 u(r)^2*r^2*dr =E2,3.",50,90);

  text("Minimise Total Energy:",50,110);

   text("integral (0.5*(du/dr)^2 - K*u^2 + P*u^2)*r^2*dr,",50,130);

  text("P(r)=0.5*E1red*integral_1 u(s)^2*min(1/r,1/s)*s^2*ds ",50,150);

  text("+ 0.5*integral_2 u(s)^2*(1/r)*s^2*ds for r in shell 1.",50,170);

  text("P(r) similar in shell 2,3.",50,190);

  text("Time step du/dt = 0.5*(d2u/dr2+2*r*du/dr)-K(r)*u+2*P*u in shell1,2,3", 10,220);

  text("+ homogeneous Neumann for shell1,2,3 at free boundary + normalisation.", 10,240);

  text("Free boundary between shell 1,2,3 established by continuity of u(r).",10,260);

  text("red: computed solution u1(r)",40,280);

  text("energy shell1 = ",40,320);

  text("energy shell2 = ",40,335);

  text("energy shell3 = ",40,350);

  text("total energy = ",40,365);

  text("r-axis",350,210);

  text("u-axis",10,10);

  text("free boundary1 =",40,380);

  text("picometer",230,380);

  text("au",180,380);

  text("picometer",230,395);

  text("au",180,395)

  text(floor(1000*M1*h)/1000,150,380);

  text(floor(1000*M1*h*0.529)/1000,200,380);

  text("free boundary2 =",40,395);

  text(floor(1000*M2*h)/1000,150,395);

   text(floor(1000*M2*h*0.529)/1000,200,395);

  text("normu1 =",40,300);

  text("normu2 =",230,300);

  text("normu3 =",230,280);

  ellipse(0,200,10);

  text(floor(10000*normu1)/10000,100,300);

  text(floor(10000*normu2)/10000,290,300);

  text(floor(10000*normu3)/10000,290,280);

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

  fill(255,0,0,255);

  ellipse(i,200-2*pow(u[i],2)*pow(i*h,2),4);

  fill(0);

  ellipse(i,200,2);

  }

  text(floor(10000*energy1)/10000,150,320);

  text(floor(10000*energy2)/10000,150,335);

  text(floor(10000*energy3)/10000,150,350);

  text(floor(10000*(energy1+energy2+energy3))/10000,150,365);

stroke(0);

line(5,0,5,200);

}

}

function startsimulation(){

  S=1;

}

function stopsimulation(){

  S=0;

}