//mode solver with sine series expansion
//yves moreau, may 98
import java.awt.*;

class   s_solver extends Object{   
// protected:
    modAppl appl;
    drawing pane;    
    rect dom;  //  domain & beginning of list of shapes
    int    mx, mxe; // sine numbers in x
    int    my, mye; // sine numbers in y
    double wl;      // wavelength
    int    nmodes;  // # of modes found
    int    mo;      // mode offset
    int    nm;      // matrix order
    double m[][];     // system matrix
    double ne[];     // effective indices
    double kx[];     // Spatial frequencies in x
    double ky[];     // Spatial frequencies in y

public s_solver(modAppl app, int mxe, int mxo, int mye, int myo)
{ this.appl = app;
  this.pane = app.pane;
  this.mx = mxe + mxo;
  this.mxe = mxe;
  this.my = mye + myo;
  this.mye = mye;
  this.wl = this.pane.lambda;
  this.nmodes = 0;  
  this.mo = 0;//mode offset
  
  this.nm = mx * my;
  try { m = new double[nm+1][nm+1];} catch (Exception e) {}; //dmatrix(1, nm, 1, nm);
  for ( int i = 1; i <= nm; i++ )
    { m[i][0] = i; for ( int j = 1; j <= nm; j++ )
      m[i][j] = 0.0;
     }
  try {ne = new double[1+nm]; } catch  (Exception e) {}; //dvector(1, nm);  
  try { kx    = new double[nm+1];} catch (Exception e) {}; // assert( kx != nil );
  try { ky    = new double[nm+1];} catch (Exception e) {}; // assert( ky != nil );
  this.dom =  new rect(appl.last_(),pane);
  this.dom.setRect(pane.centerX - pane.halfX, pane.centerX +pane.halfX,pane.centerY- pane.halfY,pane.centerY+ pane.halfY,pane.index);
  BuildSpatialFrequencies();
    } // end constructor
public rect domain_() { return(this.dom);}
private void BuildSpatialFrequencies( )  {
      double pix = Math.PI / dom.rect().width;
      double piy = Math.PI /  dom.rect().height;
      for ( int i = 0 ; i < mx ; i++ )
            for ( int j = 0 ; j < my ; j++ ) {
                  int k = my * i + j;
                  kx[k] = pix * ( i < mxe ? 2*i+1 : 2*(i+1-mxe) );   
                  ky[k] = piy * ( j < mye ? 2*j+1 : 2*(j+1-mye) );
                  }
    }
private void BuildLinearSystem( )    {
  double Lx =dom.rect().width; // domain width
  double Ly =dom.rect().height; // domain height
  double kf = 4.0 / Lx / Ly;
  // Loop over all elements of the upper-half of the system matrix 'm'
  for ( int i = 1; i <= nm ; i++ ) {
    for ( int j = 1 ; j <= i ; j++ ) {
      shape tt = dom.prev_(); while (tt instanceof shape) {  
             double ss = tt.s_integral(kx[i-1],kx[j-1],ky[i-1],ky[j-1]);
            if (tt==dom) appl.displayMsg("\ni :"+i+",j:"+j+",kx1:"+kx[i-1]+",kx2:"+kx[j-1]+",ky1:"+ky[i-1]+",ky2:"+ky[j-1]+", ss:"+ss);
            m[i][j] += kf *tt.index  * ss;// ssx * ssy ;
            tt = tt.prev_();     
      }      
      // This matrix is symmetric...
      m[j][i] = m[i][j];    
    }
  }

  double k0 =wl/2.0/Math.PI; 
  k0 *= k0; 
  for ( int i = 1 ; i <= nm ; i++ )
    {m[i][i]   -= k0 * ((kx[i-1]*kx[i-1])*pane.scaleX*pane.scaleX +(ky[i-1]*ky[i-1])*pane.scaleY*pane.scaleY );
     m[i][i]   += kf*dom.index*dom.Area()/4;
    }
}

void unPrepareShapes( )
{ shape tt=dom.prev_();
  while (tt instanceof shape) { 
    tt.index = Math.sqrt(tt.index+ dom.index); tt = tt.prev_();}
    dom.index = Math.sqrt(dom.index);
}
public int solve( boolean dispMat )
{ double nmin = dom.index; double nmax=nmin;
	dom.index = (dom.index*dom.index);
	shape tt=dom.prev_();
  while (tt instanceof shape) { 
  	if (tt.index < nmin) nmin= tt.index;
  	if (tt.index > nmax) nmax= tt.index;
    tt.index = (tt.index*tt.index) - dom.index; tt = tt.prev_();} 
 if (appl.nmin !=0) nmin= appl.nmin;
 if (appl.nmax !=0) nmax= appl.nmax;
 appl.displayMsg("\nnmin"+nmin+",nmax "+nmax);
 if (nmin >= nmax) return(nmodes);
 double  e[] = new double[1+ nm];
  Point p = dom.basePoint();
  tt=this.dom; while (tt instanceof shape) {  
  tt.offset( -p.x, -p.y ); tt= tt.prev_();}
 BuildLinearSystem(); appl.displayMsg("\nlinear system ("+nm+"x"+nm+") built,...");
 if (dispMat) appl.displayMat(m,nm);
 matlab.tred2( m, nm, ne, e );    // Householder reduction
 appl.displayMsg("\nHouseholder reduction done, ...");
 matlab.tqli( ne, e, nm, m );     // eigenvalue calculation
 appl.displayMsg("\n"+nm+" eigen values computed...");
 matlab.eigsrt(ne, m, nm);        // eigenvalue sort (descending)
 appl.displayMsg(", and sorted"); 
  mo = -1;
  for ( int i = 1; i <= nm; i++ )
    { if ( ne[i] >(nmin*nmin) && ne[i] <(nmax*nmax) )
             {nmodes++;  if ( mo == -1 )  mo = i;}
      }
  tt=this.dom; while (tt instanceof shape) {  
  tt.offset( p.x, p.y ); tt= tt.prev_();}
  unPrepareShapes();
  e=null; //  free_dvector(e, 1, nm);  
  return nmodes;
}

public s_Mode getMode( int n )
    {double e[];
      if ( nmodes == 0 ) return null;
      if ( n >= nmodes ) return null;
      n += mo;
     try {
       e = new double[nm];
     	} catch (Exception ex) { return null;};
    for ( int i = 1; i <= nm; i++ ) e[i-1] = m[i][n];
      return new s_Mode( this, Math.sqrt(ne[n]), wl, e );
    }
public double EffectiveIndex( int n )
    {  if ((n+mo)<0) return(0);
      return Math.sqrt(ne[n + mo]); }
public void finalize() { m=null; ne=null; kx=ky=null;}
}// end class s_solver