// This is the code which is associated with the web page
//   http://www.fen.bilkent.edu.tr/~yalabik/applets/1d_bound.html
// Original code was in Java (circa 2000), this is a simpler version in Javascript
// Cemal Yalabık - 2025  yalabik@fen.bilkent.edu.tr 


  const canvas = document.getElementById('drawCanvas');
  const ctx = canvas.getContext('2d');
  
  const neg_v = canvas.height;
  const pos_v = canvas.height * 2 / 3;
  const def_v = pos_v + canvas.height / 6;
  const psi_zero = canvas.height / 3;
  const p_base = canvas.height/3;
  
  const n = canvas.width;
  const m = 5 ; // no of eigenvectors
  const scale = 0.0001 ; // plot to potential scale
  const p_scale = 0.1*(neg_v-pos_v)*Math.sqrt(n);  // scales avefunction plot
  const pot_shift = 0. // 1.; // eigenvalue shift in Hamiltonian
  
  var pot  = Array(n).fill(0);
  var diag = Array(n).fill(0);
  var rhs  = Array(n).fill(0);
  var prev_rhs = Array(n).fill(0);
  var eig  = Array(m).fill(0);
  
  var xbef = 0;
  var old_pot=0;
  
  default_pot();
  
  solve();

  let drawing = false;
  var solving = false;
  var pot_changed = true;

  canvas.addEventListener('mousedown', e => {  //////////////////////
    if((! solving) && -1 < e.offsetX && e.offsetX < n && e.offsetY < neg_v && e.offsetY > pos_v){
      drawing = true;
      pot[e.offsetX] = e.offsetY;
      xbef = e.offsetX;
      draw_pot();
      pot_changed = true;
    }
     
  });

  canvas.addEventListener('mousemove', e => {  //////////////////////

//    ctx.font = "20px Arial";
//    ctx.clearRect(0, 0, canvas.width, 100);

    if((! solving) && drawing && -1 < e.offsetX && e.offsetX < n && e.offsetY < neg_v && e.offsetY > pos_v){
      old_pot = pot[xbef];
//      ctx.fillText("xbef = " + xbef + "  pot[xbef] = " + pot[xbef],200,50);
      if(xbef < e.offsetX){
        for( let i=xbef;i<= e.offsetX;i++ ){
          pot[i] = old_pot + (i-xbef)*(e.offsetY-old_pot)/(e.offsetX-xbef);
//          ctx.fillText("xbef = " + xbef + " old_pot = " + old_pot + " e.offsetX = " + e.offsetX + " e.offsetY = " + e.offsetY ,10,20);
//          ctx.fillText("i = " + i + " pot[i] = " + pot[i],10,50);
        }
      }
      else if(xbef > e.offsetX){
        for( let i=xbef;i>= e.offsetX;i-- ){
          pot[i] = old_pot + (i-xbef)*(e.offsetY-old_pot)/(e.offsetX-xbef);
        }        
      }
      xbef = e.offsetX;
//      ctx.fillText("out!",500,50);
      draw_pot();
      pot_changed = true;
    }
    else{ drawing = false;}
  });

  canvas.addEventListener('mouseup', () => {  //////////////////
    drawing = false;
    if( ! solving ) solve();
  });

  canvas.addEventListener('mouseleave', () => {   ///////////////////
    drawing = false;
    if( ! solving ) solve();
  });

  function clearCanvas() {   //////////////////////
    ctx.clearRect(0, 0, canvas.width, canvas.height);
  }
  
  function default_pot(){  ////////////////
    for(let i=0;i<n;i++){
      if(i<n*9/20 || i> n*11/20) pot[i] = (pos_v+def_v)*0.5;
      else pot[i] = neg_v -5;
    }
    draw_pot();
    pot_changed = true;
  }
  
  function draw_pot(){  //////////////////////////////
      ctx.clearRect(0, pos_v, n, neg_v);
      ctx.strokeStyle = 'black';
      ctx.beginPath();  
      ctx.moveTo(0, pot[0]);    
      for (let i = 1; i < n; i++) {
        ctx.lineTo(i, pot[i]);
      }
      ctx.stroke();
  }
  
  function solve(){  ////////////////////////////////
  
// n points, m eigenvectors 
// We are interested in eigenvalues and eigenvectors of -D^2 + V
// Smallest eigenvalues are largest eigenvalues of 1/(-D^2 + V) 
// We do repeated multiplications U(k+1) = 1/(-D^2 + V) U(k)
// which is equivalent to solving (-D^2 + V) U(k+1) = U(k)
// Matrix elements look like 
// ...
//    -1   2+V  -1
//         -1   2+V  -1
//              -1   2+V  -1
//                          ...
// A constant (pot_shift) may be needed to be added to the diagonal
//  to give a positive shift if lowest eigenvalue may go negative.
  
  if( ! pot_changed ) return;
  
  solving = true;
  
  let m_conv = 20; // min no of loops 
  let size = 0.;
  let xsize = 0.;
  let ssize = 0.;
  let yy = 0.;
  let error = 0.;
  
  let rhs = Array(n).fill(0);
  let prev_rhs = Array(n).fill(0);
  let rhso = new Array(m); // Creates a new array with m undefined elements
    for (let i = 0; i < m; i++) {
      rhso[i] = new Array(n).fill(0); // Each sub-array has n elements filled with 0
    }
    
  
  for(let mloop=0;mloop<m;mloop++){
 
      // fill in initial random vector
    for(let i=0;i<n;i++){ 
      rhs[i]=Math.random();
      prev_rhs[i]=0;
    }
      
    for(let loop=0;loop<m_conv;loop++){

      for(let j=0;j<mloop;j++){  // project out previous solutions
        let proj= 0.;
        for(let i=0;i<n;i++){
          proj+=rhso[j][i]*rhs[i]; // the inner product (normalized solutions are in rhso)
        }
        for(i=0;i<n;i++){
          rhs[i]-=proj*rhso[j][i];  // subtract it out
        }
      }

// working on the mloop'th vector

      size = 0.;  // size of vector before operation
      for(let i=0;i<n;i++){
        size+=rhs[i]*rhs[i];
      }


      for(let i=0;i<n;i++){
        diag[i] = pot_shift + scale*(neg_v-pot[i]) + 2.;  // The diagonal elements
      }

        
// the tridiagonal equation solution:
      for(let i=0;i<n-1;i++){
        let q=-1./diag[i];
        diag[i+1]=diag[i+1]+q;
        rhs[i+1]=rhs[i+1]-rhs[i]*q;
      }
      rhs[n-1]=rhs[n-1]/diag[n-1];
      for(let i=n-2;i>=0;i--){
       rhs[i]=(rhs[i]+rhs[i+1])/diag[i];
      }
// end of the tridiagonal equation solution

// normalize the vector
      xsize=0;  // size of vector after operation
      for(let i=0;i<n;i++){
        xsize+=rhs[i]*rhs[i];
      }
      ssize= 1./Math.sqrt( xsize) ;
      biggest=0.;
      for(let i=0;i<n;i++){
        rhs[i]*=ssize;    //  rhs normalized
        if(Math.abs(rhs[i]) > biggest){ biggest=Math.abs(rhs[i]); }
      }
      if(biggest != 0.){ biggest=1./biggest;}  // contains scale max to 1
      biggest = 1. // testing no scaling
      
      eig[mloop]= Math.sqrt( size/xsize )-pot_shift; // present estimate

// plot curves at each iteration
      ctx.clearRect(0, 0, n, pos_v);
      for(let i=0;i<=mloop;i++){
        if(i==0){ ctx.strokeStyle = 'red';}
        if(i==1){ ctx.strokeStyle = 'blue';}
        if(i==2){ ctx.strokeStyle = 'green';}
        if(i> 2){ ctx.strokeStyle = 'black';}
        
        if(i<mloop && boxes[i] == 1){
          ctx.beginPath();  
          ctx.moveTo(0, p_base-p_scale*rhso[i][0]);    
          for (let j = 1; j < n; j++) {
            ctx.lineTo(j, p_base-p_scale*rhso[i][j]);
          }
        }else if(boxes[mloop] == 1){
          ctx.beginPath();  
          ctx.moveTo(0, p_base-p_scale*rhs[0]*biggest);    
          for (let j = 1; j < n; j++) {
            ctx.lineTo(j, p_base-p_scale*rhs[j]*biggest);
          }
        }
        
        ctx.stroke();
      }
      

      error=0.
      for(let i=1;i<n;i++){
        error+=Math.abs(rhs[i]-prev_rhs[i]);
        prev_rhs[i]=rhs[i];
      }

      if(error*p_scale < 0.1 ) { break;}

    }  // ends loop

// mloop'th eigenvalue and eigenvector determined.
// Plot eigenvalue level  
    if(mloop==0){ ctx.strokeStyle = 'red';}
    if(mloop==1){ ctx.strokeStyle = 'blue';}
    if(mloop==2){ ctx.strokeStyle = 'green';}
    if(mloop> 2){ ctx.strokeStyle = 'black';} 

    yy= neg_v - eig[mloop]/scale;
    if(yy>pos_v && yy<neg_v){
      ctx.fillText("eigenvalue " + mloop + " = " + (eig[mloop]/scale).toPrecision(4),3*n/4+5,yy);
      ctx.moveTo(n/4, yy);    
      ctx.lineTo(3*n/4, yy);
      ctx.stroke();
    }

    for(let i=0;i<n;i++){
      rhso[mloop][i]=rhs[i];
    }
    
  }  // closes mloop

  pot_changed = false;  
  solving = false;
  
 }
 
 function checkbox_click(){   ////////////////////////
   var box = document.getElementById("lev0");
   if (box.checked == true) boxes[0]=1;
   else boxes[0]=0;
   box = document.getElementById("lev1");
   if (box.checked == true) boxes[1]=1;
   else boxes[1]=0;
   box = document.getElementById("lev2");
   if (box.checked == true) boxes[2]=1;
   else boxes[2]=0;
   box = document.getElementById("lev3");
   if (box.checked == true) boxes[3]=1;
   else boxes[3]=0;
   box = document.getElementById("lev4");
   if (box.checked == true) boxes[4]=1;
   else boxes[4]=0;
   
   if( ! solving ){pot_changed = true; solve();}
   
 }