/**
 *
 * Francois Nedelec, 2001
 * EMBL, Heidelberg, Europe
 * nedelec@embl-heidelberg.de
 * 
 */

import java.awt.*;
import java.awt.event.*;
import java.util.Random;
import java.text.NumberFormat;


class Aster extends Canvas implements Runnable  {
    //create a square canvas with an aster in the center;
    
    //the thread of the simulation, used to start and stop it
    Thread simulatorThread = null;
    
    //the size of the canvas area
    public static final int wSize = 600;
    
    //a graphic buffer used for drawing smoothly
    Image bi = null;
    Graphics big;

    //used to format numbers:
    NumberFormat nf = NumberFormat.getInstance();
 
    //a random number generator:
    Random RNG = new Random();

    //parameters of the simulation:
    int model=0, nbMTs = 0, acceleration = 1, clock = 10 ;
    long time;
    double zoom;       // zoom in space
    double zoomN = 1;  // zoom for the y-axis of the length-distribution
    double dt;         // time step in seconds
    double timesim;    // simulated time so far, seconds;
    double meanlength, theorylength; // 

    
    //parameters of the model
    double growSp, cataFq, shrinkSp, rescFq;
    double grow, shrink, cata, resc;    //grow = dt * growSp,  etc.
    
    //state variables of the microtubules:
    double[] MTlength;
    int[] MTstate;

    //we precompute the microtubule's directions (faster drawing)
    double[] MTcos;
    double[] MTsin;

    //variables for the curve of length distribution:
    final int maxprofile = 60;       // size of profile arrays
    final int sumcnt = 138;          // number of profiles recorded
    final int binL = 2;              // each profile point is 1/binL micro-m

    double avginterval = 100;        // time interval used to compute profile
    
    // the profile profile[ x ] counts how many microtubules are 
    // between length x / binL and ( x+1 ) / binL.
    // the second index is the profile number, related to time
    long[][] profile = new long[maxprofile][sumcnt];
    
    // the cumulative sum of the profiles
    long[] sumprofile = new long[maxprofile];
    
    // pindx is the current profile number
    int binT = 1, pcnt = 0, pindx = 0;
    
    // used in drawing the profile
    int[] xPoints = new int[ maxprofile ];
    int[] yPoints = new int[ maxprofile ];
    


    public Aster() {
        //global initialization:
        nf.setMaximumFractionDigits(0);
        nf.setMinimumFractionDigits(0);
        setSize(wSize, wSize);
    }
   
        
    //set the number of filaments
    public void setNumber( int mt ) {
        
        if ( ( mt > 0 ) && ( mt != nbMTs ) ) {
            //reset the arrays if we change the number of filaments:
            MTlength = new double[mt];
            MTcos    = new double[mt];
            MTsin    = new double[mt];
            MTstate  = new int[mt];

            //recalculate and store the angles
            double angle;
            for(int i=0; i < mt; ++i) {
                angle = 2 * i * Math.PI / mt;
                MTcos[ i ] = Math.cos( angle ); 
                MTsin[ i ] = Math.sin( angle );
            }
            nbMTs = mt;
            reset( 0 );
        }
    }
    
    public void reset(int s) {
        timesim = 0;
        //reset the microtubule with state s given as argument:
        for( int i=0; i < nbMTs; ++i) {
             MTstate[i]  = s; 
             MTlength[i] = 0; 
        }
        // reset the stored profiles:
        pindx = 0;
        for( int i=0; i < maxprofile; ++i ) {
           sumprofile[i] = 0;
           for( int p=0; p < sumcnt; ++p ) 
                profile[i][p] = 0;
        }
    }
    	
    ///set the zoom for drawing in X and Y
    public void setZoom( int Z ) {
        zoom  = Z / 10.0;
    }
    
    ///set the zoom for the MT count profile
    public void setZoomN( int Z ) {
        zoomN = Z / 1000.0;
    }
    
    public void setTiming(int accl, int clockrate, int avg ) {
      timesim       = 0;
      acceleration  = accl;
      avginterval   = avg;
      clock         = clockrate;
      precomputeParams();
    }
        
    private void precomputeParams() {
        dt = 0.001 * acceleration * clock;
      //System.out.println( nf.format( dt ) );
        grow   = dt * growSp;
        shrink = dt * shrinkSp;
        cata   = dt * cataFq;
        resc   = dt * rescFq;
        if ( ( cata > 0.2 ) || ( resc > 0.2 ) ) {   
          //JOptionPane.showMessageDialog(this, "acceleration too big");
        }
        
        // calculate binT, the frequency at which profiles are stored:
        binT = (int)Math.round( avginterval / ( sumcnt * dt ) );
        if ( binT < 1 ) binT = 1;
        //System.out.println( nf.format( binT ) );
    }

    public void setParams(double G, double S, double C, double R) {
        //rates/speeds are given per minutes, we transform to per seconds
        rescFq   = R / 60.0;
        shrinkSp = S / 60.0;
        cataFq   = C / 60.0;
        growSp   = G / 60.0;
        //predicted length for the 'classic' model:
        theorylength = growSp * shrinkSp / ( shrinkSp * cataFq + growSp * rescFq );
        precomputeParams();
    }
    
    ///set the type of model for MT dynamic instability
    public void setModel( int m ) {
        model = m;
    }
      
    // main graphic procedure
    public void paint(Graphics g) {
        
      //set up the graphical buffer, on which the image is drawn first
      if ( bi == null ) {
        bi  = createImage(wSize, wSize);
        big = bi.getGraphics();
      }
      
      //======================paint surface in black:
      big.setColor(Color.black);
      big.fillRect(0, 0, wSize, wSize);
      
      //======================little clock:
      final int clocksize = 30;
      int x = wSize - clocksize - 5;
      int y = clocksize + 20 ;
      
      big.setColor(Color.yellow); 
      big.drawArc( x-clocksize, y-clocksize, 2*clocksize, 2*clocksize, 0, 360 );
      
      //yellow hand for real time: 
      double d = System.currentTimeMillis() * acceleration * Math.PI / 30000.0;
      big.drawLine( x, y, x+(int)Math.round(clocksize*Math.cos(d)),	
                    y+(int)Math.round(clocksize*Math.sin(d)));
      
      //white hand for simulation time:
      big.setColor(Color.white);
      d = time * acceleration * Math.PI / 30000.0;
      big.drawLine( x, y, x+(int)Math.round(clocksize*Math.cos(d)),
                    y+(int)Math.round(clocksize*Math.sin(d)));
      
      nf.setMaximumFractionDigits(0);
      big.drawString( nf.format( timesim )+" s", x - clocksize, 15 );
      

      //================================scale bar:
      final int left = 20;
      
      big.setColor(Color.white);
      y = (int)Math.min( 5, zoom );
      big.fillRect(left, 5, (int)Math.round( 10 * zoom ), y);
      x = left + (int)Math.round( 11 * zoom );
      big.setColor(Color.blue);
      big.drawString("micro-meters", x, y+7 );
      
      //=================little ticks every micro-meter
      big.setColor(Color.black);
      for(int i=1; i<10; ++i) {
        x = left + (int)Math.round( i * zoom );
        big.drawLine( x, 5, x, 5+y );
      }
      
      //================small vertical bar at mean length
      big.setColor(Color.white);
      x = left + (int)Math.round(zoom * meanlength);	
      big.drawLine(x, 12, x, 18);
      
      //================small vertical bar at the theoretical mean length:
      big.setColor(Color.yellow);
      if ( theorylength > 0 ) {
        x = left + (int)Math.round(zoom * theorylength);
        big.drawLine(x, 12, x, 18);
        nf.setMinimumFractionDigits(2);
        big.drawString(nf.format(theorylength), x+5, 19 );
      }
      else
        big.drawString("unbounded growth", left, 19 );

      int center_x = wSize / 2;
      int center_y = wSize / 2;
      
      //draw the microtubules:
      for(int i=0; i < nbMTs; ++i ) 
        if ( MTlength[i] > 0 ) {
          d = zoom * MTlength[ i ];
          y = center_y + (int)Math.round( MTsin[i] * d );
          x = center_x + (int)Math.round( MTcos[i] * d );
          if ( MTstate[i] == 0 ) 
            big.setColor(Color.white);
          else
            big.setColor(Color.yellow);
          big.drawLine(center_x, center_y, x, y );
        }
          

      //========================the profile of length distribution:
        
      double ZL  = zoom / binL;
      double ZN  = zoomN * binL * ( wSize - 50 );
      double ZNS = ZN / sumcnt;
      int dx     = (int) Math.round( ZL );
      int maxp   = (int) Math.min( maxprofile-1, Math.round((wSize-50)/ZL));
      long maxval;
      
      //draw the axis:
      big.setColor(Color.blue);
      big.drawLine(left, 19, left+(int)Math.round( maxp*ZL ), 19);
      big.drawLine(left, 19, left, wSize-25);
      //y-axis labels:
      
      for( int i=1; i < 10 ; ++i ) {
        y = 20+(int)Math.round( ZN * i ) ;
        big.drawLine(left-10 , y, left, y );
      }
      big.drawString(nf.format( 5 ), 1 , 20+(int)Math.round( ZN * 5 ) );
      for( int i=10; i < Math.min( nbMTs, 100) ; i+= 10 )
        big.drawString(nf.format( i ), 1 , 20+(int)Math.round( ZN * i ) );
      big.drawString("MT / um", 1 , wSize - 10 );  
      
      //draw the curve for instantaneous MT count
      big.setColor(Color.yellow);	
      for(int i=0; i < maxp; ++i) {
        x = left + (int) Math.round( ZL * i );
        xPoints[ i ] = x;
        y = 20 + (int) Math.round( ZN * profile[ i ][ pindx ] );
        yPoints[ i ] = 20 + (int) Math.round( ZNS * sumprofile[ i ] );
        if ( y < wSize )  big.drawLine( x, y, x+dx, y );
      }

      //draw the curve for averaged MT count
      big.setColor(Color.white);
      big.drawPolyline(xPoints, yPoints, maxp );
      
      //copy the buffered image on the screen
      g.drawImage(bi, 0, 0, Color.black, null);
    }

    
    public void update(Graphics g) {
      paint( g );
    }
    

    //classical model: MT transitions are independent of length
    private void iterateClassic() {
        
      for(int i=0; i < nbMTs; ++i) {
        switch( MTstate[i] ) {
          case 0:
            MTlength[i] += grow;
            if ( RNG.nextDouble() < cata ) 
              MTstate[ i ] = 1;
              break;
          case 1:
            MTlength[i] += shrink;
            if ( RNG.nextDouble() < resc ) 
              MTstate[ i ] = 0;
              break;
        }
      }
    }
    
    //alternate phenomenological model: transition dependent on the length of the MT
    private void iterateModern() {
        
      double cataM = cata / 10.0;
      double rescM = resc / 3.0;
      
      for(int i=0; i < nbMTs; ++i) {
        switch( MTstate[i] ) {
          case 0:
            MTlength[i] += grow;
            if ( RNG.nextDouble() < cataM * MTlength[i] )
              MTstate[ i ] = 1;
              break;
          case 1:
            MTlength[i] += shrink;
            if ( RNG.nextDouble() < rescM * ( 13 - MTlength[i] ) ) 
              MTstate[ i ] = 0;
              break;
        }
      }
    }
    
    
    ///perform one monte-carlo iteration
    private void iterate() {
	
      time += clock;
      timesim += dt;
      meanlength = 0;
      int p;
	
      //calculate the sum of all profiles:
      if ( ++pcnt >= binT ) {
        pcnt  = 0;
        int pnext = ( pindx + 1 ) % sumcnt;
        for( p = 0; p < maxprofile; ++p )
          sumprofile[p] += profile[p][ pindx ] - profile[p][ pnext ];
        pindx = pnext;
        //if (pnext == 0) System.out.println("time "+nf.format(timesim));
      }
      
      for( p=0; p < maxprofile; ++p )
        profile[ p ][ pindx ] = 0;
      
      if ( model == 1 ) 
        iterateModern();
      else
        iterateClassic();
      
      //calculate the mean-length, and set the length distribution
      for(int i=0; i < nbMTs; ++i) {
        
        if ( MTlength[i] <= 0 ) {  
          MTlength[i] = 0;
          MTstate[i]  = 0;        //rescure MT of null length
        }
        
        meanlength += MTlength[ i ];
        
        p = (int)Math.floor( binL * MTlength[i] );
        if ( p >= maxprofile ) p = maxprofile - 1;
        ++profile[ p ][ pindx ];
      }
      
      meanlength /= nbMTs;
    }
    
    ///run the simulator thread
    public void run() {
      timesim = 0;
      time = System.currentTimeMillis();
      while( Thread.currentThread() == simulatorThread ) {
        iterate();   
        repaint();
        try {
          long delay = Math.max( 0, time - System.currentTimeMillis() ); 
          simulatorThread.sleep( delay );      
        }
        catch( InterruptedException e ) { return; }
      }
    }   
    
    ///create a new thread
    public void start() {
        //System.out.println("aster starts");
      if ( simulatorThread == null ) 
        simulatorThread = new Thread( this );
      simulatorThread.start();
    }
    
    ///stop the thread
    public void stop() {
      simulatorThread = null;
    }
    
}