createCanvas(400, 400);

}

​

function draw() {

  background(220);

}

​

let u = [];

let K = [];

let Z;

let E1=2;

let E1red=1/2;

let E2;

let E3;

let E2red;

let E3red;

let N=400;

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 start;

let P = [];

let fix=0;

let D=5;

let S=0;

let diff=0.001;

​

function setup() {

  createCanvas(400, 400);

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

  slider2=createSlider(0,4,3,1);

  slider3=createSlider(0,4,3,1);

  start=createButton('Start');

  start.position(350,370);

  start.mousePressed(startsimulation);

  h=D/N;

  Z=slider1.value();

  E2=slider2.value();

  E3=slider3.value();

  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=1;i<N+1;i++){

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

  }

  for (var i=1; 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=slider1.value();

  E2=slider2.value();

  E3=slider3.value();

  noStroke();

  fill(0);

  text("2",205,200);

  text("E2",240,200);

  text("E3", 290,200);

  text("RealQM Atom Simulator: 3 Shell Structure",100,20);

  text("Kernel Charge 2 < Z < 11.",100,40);

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

  text("Choose Z, E2 and E3 and press Start",100,90);

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

  text("Lithium Z=3: 2+1",50,270);

  text("Beryllium Z=4: 2+2",50,285);

  text("Boron Z=5: 2+2+1",50,300);

  text("Carbon Z=6: 2+2+2",50,315);

  text("Nitrogen Z=7: 2+3+2",50,330);

  text("Oxygen Z=8: 2+3+3",50,345);

  text("Fluorine Z=9: 2+4+3",50,360);

  text("Neon Z=10: 2+4+4",50,375);

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

  text(Z,70,395);

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

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

  text(E2,210,395);

  text(E3,310,395);

  if (E2==0){

  M1=400;

  }

  if (E3==0){

  M2=400;

  }

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

  text("(Li observed= -7.48)",285,360);

  fix=0.1;

  D=5;

  }

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

  text("(Li+ observed= -7.28)",285,360);

  fix=0.2;

  D=5;

  }

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

  text("(Be observed= -14.47)",280,360);

  fix=0.1;

  D=5;

  }

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

  text("(Be+ observed= -13.65)",280,360);

  fix=0.4;

  D=5;

  }

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

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

  fix=0.1

  D=7;

  }

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

  text("(Bo2+ observed= -24.53)",285,365);

  fix=0.7;

  D=5;

  }

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

  text("(C observed= -37.7)",285,360);

  fix=0;

  D=7;

​

  }

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

  text("(N observed= -54.6)",285,360);

  fix=0;

  D=7;

  }

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

  text("(N+ observed= -54.0)",285,360);

  fix=0.60;

  D=7;

  }

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

  text("(observed = -74.8)",285,360);

  fix=0;

  D=7;

  }

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

  text("(O6+ observed = -59.2)",280,360);

  fix=0.8;

  D=7;

  }

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

  text("(Fluor = -99.8)",285,360);

  fix=0;

  D=7;

  }

  if (Z==9 && E2==0 && E3==0){

  text("(Fluor7+ = -75.5)",285,360);

  fix=1;

  D=7;

  }

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

  text("(Neon observed = -128.5)",285,360);

  fix=0;

  D=7;

  }

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

  text("(Neon8+ = -94.0)",285,360);

  fix=1.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 observed= -7.48)",285,360);

//  fix=0;

//  D=5;

  }

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

  text("(Li+ observed= -7.28)",285,360);

  fix=0.2;

  D=5;

  }

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

  text("(Be observed= -14.47)",280,360);

//  fix=0.4;

//  D=5;

  }

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

  text("(Be+ observed= -13.65)",280,360);

  fix=0.4;

  D=5;

  }

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

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

  fix=0.3

  D=7;

  }

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

  text("(Bo2+ observed= -24.53)",285,365);

  fix=0.7;

  D=5;

  }

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

  text("(C observed= -37.7)",285,360);

  fix=0.4;

  D=7;

​

  }

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

  text("(N observed= -54.6)",285,360);

  fix=0.55;

  D=7;

  }

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

  text("(N+ observed= -54.0)",285,360);

  fix=0.60;

  D=7;

  }

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

  text("(O observed = -74.8)",285,360);

  }

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

  text("(O6+ observed = -59.2)",280,360);

  fix=0.8;

  D=7;

  }

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

  text("(Fluor = -99.8)",285,360);

  }

  if (Z==9 && E2==0 && E3==0){

  text("(Fluor7+ = -75.5)",285,360);

  fix=1;

  D=7;

  }

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

  text("(observed = -128.5)",285,360);

  }

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

  text("(Neon8+ = -94.0)",285,360);

  }

  noStroke();

  fill(0);

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

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

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

  text("Normalisation: 0<E2<5, 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(M1*h,150,380);

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

  text(M2*h,150,395);

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

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

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

  ellipse(0,200,10);

  text(normu1,100,300);

  text(normu2,290,300);

  text(normu3,290,280);

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

  fill(255,0,0,255);

  ellipse(i,200-50*u[i],4);

  fill(0);

  ellipse(i,200,2);

  }

  text(energy1,150,320);

  text(energy2,150,335);

  text(energy3,150,350);

  text(energy1+energy2+energy3,150,365);

stroke(0);

line(5,0,5,200);

}

}

function startsimulation(){

  S=1;

}