/** ‚QŽŸŒ³”g“®•û’öŽ®�i‰~—̈æ�E”¼‰AƒXƒL�[ƒ€�j
 *
 *  <applet code = Wave2d_i2 width=500 height=500></applet>
 *
 */

import java.applet.*;
import java.awt.*;
import java.awt.event.*;
import corejava.*;
import Mitsui.*; 
import cern.jet.math.*;  

public class Wave2d_i2 extends Applet{
    int nfunc=0;
    MitsuiWorld_2 m = new MitsuiWorld_2();
    double xmin,xmax,ymin,ymax,zmin,zmax;
    private TextArea ta;

    public void init(){
        Graphics g=getGraphics();
        m.setGraphics(g);
        m.setScreenSize(getSize().width,getSize().height);
        //[a,b]�~[c,d]�~[e,f]
        xmin=-2.0;
        xmax=2.5;
        ymin=-2.0;
        ymax=2.5;
        zmin=-1.0;
        zmax=1.0;
        m.setArea2(xmin,xmax,ymin,ymax,zmin,zmax);
        m.setColor(3);
        add(ta = new TextArea(5,30));
        //ta.setBounds(0,0,10,7);
    }
    
    public void paint(Graphics g){
        int nr = 20;
        int nt = 100;
        double tau = 0.045;
        double Tmax = 3.0;
        double dt = 0.05;
        double hr = 1.0/nr;
        double ht = 2.0*Math.PI/nt;
        double lambda_r = tau / hr;
        double lambda_t = tau / ht;
        double lambdar2 = lambda_r * lambda_r;
        double lambdat2 = lambda_t * lambda_t;
        
        double[][] u    = new double[nr+1][nt+1];
        double[][] u1   = new double[nr+1][nt+1];
        double[][] newu = new double[nr+1][nt+1];
        double[][] al   = new double[nr][nt];
        double[][] ad   = new double[nr][nt];
        double[][] au   = new double[nr][nt];
        double[][] ab   = new double[nr][nt];
        double[][] ar   = new double[nr][nt];
        double[]   b    = new double[nt+1];

        if(dt < tau){
            dt = tau;
        }
        
        double skip = Math.rint(dt/tau);

        g.clearRect(0,0,getSize().width,getSize().height);
        m.hideBirdView(u,nr,nt,1);

        //LU•ª‰ð
        for(int i=1;i<nr;i++){
            double ri = i * hr;
            double ri2= ri * ri;
            double d  = 1+2*lambdat2/ri2;
            double od = -lambdat2 / ri2;
            for(int j=0;j<nt;j++){
                ad[i][j] = d;
                al[i][j] = au[i][j] = od;
                ab[i][j] = ar[i][j] = 0.0;
            }
            ab[i][0] = ab[i][nt-2] = ar[i][0] = ar[i][nt-2] = od;
            ptrilu0(nt,al[i],ad[i],au[i],ab[i],ar[i]);
        }

        for(int i=0;i<=nr;i++){
            double ri = i*hr;
            for(int j=0;j<=nt;j++){
                double theta_j = j*ht;
                u1[i][j] = phi(ri,theta_j);
            }
        }
        for(int i=0;i<=nr;i++){
            double ri = i*hr;
            for(int j=0;j<=nt;j++){
                double theta_j = j*ht;
                u[i][j] = phi(ri,theta_j) +tau*psi(ri,theta_j);
            }
        }
        

        //ŽžŠÔƒ‹�[ƒv
        double nMax = Math.rint(Tmax / tau);
        for(int n=1;n<=nMax;n++){
            for(int i=1;i<nr;i++){
                //‰E•Ó‚ð�ì‚é
                double alpha = lambdar2 * (1.0+0.5/i);
                double beta  = lambdar2 * (1.0-0.5/i);
                for(int j=0;j<nt;j++){
                    b[j] = 2*(1 - lambdar2)*u[i][j] + alpha * u[i+1][j]
                        + beta * u[i-1][j] -u1[i][j];
                }
                //˜A—§‚PŽŸ•û’öŽ®‚ð‰ð‚­
                ptrisol0(nt,al[i],ad[i],au[i],ab[i],ar[i],b);
                //ƒRƒs�[
                for(int j=0;j<nt;j++)
                    newu[i][j] = b[j];
                newu[i][nt] = b[0];
            }
            //Dir
            for(int j=0;j<=nt;j++)
                newu[nr][j] = 0.0;
            
            //Œ´“_
            double s = 0.0;
            for(int j=0;j<nt;j++)
                s += u[1][j];
            
            newu[0][0] = 4*lambdar2 * (s/nt-u[0][0]) + 2*u[0][0] - u1[0][0];

            for(int j=1;j<=nt;j++)
                newu[0][j] = newu[0][0];
            
            // newu, u, u1 �X�V
            for(int i=0;i<=nr;i++)
                for(int j=0;j<=nt;j++)
                    u1[i][j] = u[i][j];
            for(int i=0;i<=nr;i++)
                for(int j=0;j<=nt;j++)
                    u[i][j] = newu[i][j];
            
            g.clearRect(0,0,getSize().width,getSize().height);
            m.hideBirdView(u,nr,nt,1);
                            
            double t = n * tau;
    
            //Œë�·‚𑪂é
            double maxe,e,x=0,y=0;
            maxe = 0;
            e = 0;
            for(int i=0;i<=nr;i++){
                double ri = i*hr;
                for(int j=0;j<=nt;j++){
                    double theta_j = j * ht;
                    e = Math.abs(exactu(ri,theta_j,t) - u[i][j]);
                    if (e > maxe) {
                        maxe = e;
                        x = ri * Math.cos(theta_j);
                        y = ri * Math.sin(theta_j);
                    }
                }
            }
            Format fmt = new Format("%7.2f");
            Format fmt_e = new Format("%7.2e");
            String str =ta.getText();
            if (str.length() > 1000)
                ta.setText("Œë�·"+fmt_e.form(e)+"\n(x,y)="
                           +fmt.form(x)+","+fmt.form(y)+")\n");
            else
                ta.setText(str+"Œë�·"+fmt_e.form(e)+"\n(x,y)="
                           +fmt.form(x)+","+fmt.form(y)+")\n");
            if(skip == 0){
                t = n*tau;
            }
        }
    }

    double mu01 = 2.404825557695773;
    double mu11 = 3.831705970207512;

    //�‰Šú’l
    public double phi(double r, double theta){
        return Bessel.j0(mu01*r); 
    }

    public double psi(double r,double theta){
        return Bessel.j1(mu11*r)*Math.cos(theta);
    }

    //Œµ–§‰ð
    public double exactu(double r,double theta,double t){
        return Bessel.j0(mu01*r)*Math.cos(mu01*t)+1/mu11*Bessel.j1(mu11*r)*
                Math.cos(theta)*Math.sin(mu11*t);
    } 
    
    public double maxnorm(int m,int n,double[][] u){
        int i0 = 0;
        int j0 = 0;
        double absu;
        double tmpmax = Math.abs(u[0][0]);
        for(int i=0;i<m;i++){
            for(int j=0;j<n;j++){
                if((absu = Math.abs(u[i][j]))>tmpmax){
                    tmpmax = absu;
                    i0 = i;
                    j0 = j;
                }
            }
        }
        return tmpmax;
    }

    public void ptrilu0(int n,
                        double[] al,double[] ad,double[] au,
                        double[] ab,double[] ar){
        int nm1 = n - 1;
        int nm2 = n - 2;
        for(int k=0;k<nm2;k++){
            int kp1 = k + 1;
            
            al[kp1] /= ad[k];
            ad[kp1] -= al[kp1]*au[k];
            if(kp1 < nm2)
                ar[kp1] -= al[kp1]*ar[k];
            else
                au[kp1] -= al[kp1]*ar[k];
            ab[k] = ab[k]/ad[k];
            if(kp1 < nm2)
                ab[kp1] -= ab[k]*au[k];
            else
                al[nm1] -= ab[k]*ar[k];
            ad[nm1] -= ab[k]*ar[k];
        }
        al[nm1] /= ad[nm2];
        ad[nm1] -= al[nm1]*au[nm2];
    }
    
    public void ptrisol0(int n,
                         double[] al,double[] ad,double[] au,
                         double[] ab,double[] ar,double[] b){
        int nm1 = n - 1;
        for(int k=0;k<nm1;k++){
            b[k+1] -= al[k+1]*b[k];
            if(k+1 < nm1 && k+1 > nm1)
                b[nm1] -= ab[k]*b[k];
        }
        b[nm1] = b[nm1]/ad[nm1];
        for(int k = n-2;k>=0;k--){
            b[k] -= au[k]*b[k+1];
            if(k+1 < nm1 && k+1 > nm1)
                b[k] -= ar[k]*b[nm1];
            b[k] /= ad[k];
        }
    }
    public void ptrilu1(int n,
                        double[] al,double[] ad,double[] au,
                        double[] ab,double[] ar){
        int nm1 = n - 1;
        for(int k=1;k<nm1;k++){
            int kp1 = k + 1;
            
            al[kp1] /= ad[k];
            ad[kp1] -= al[kp1] * au[k];
            if(kp1<nm1)
                ar[kp1] -= al[kp1]*ar[k];
            else 
                au[kp1] -= al[kp1]*ar[k];

            ab[k] /= ad[k];
            if(kp1<nm1)
                ab[kp1] -= ab[k]*au[k];
            else
                al[n] -= ab[k]*au[k];
            ad[n] -= ab[k]*ar[k];
        }
        al[n] /= ad[nm1];
        ad[n] -= al[n]*au[nm1];
    }

    public void ptrisol1(int n,
                         double[] al,double[] ad,double[] au,
                         double[] ab,double[] ar,double[] b){
        for(int k=1;k<n;k++){
            b[k+1] -= al[k+1]*b[k];
            if(k+1<n && k+1>n)
                b[n] -= ab[k]*b[k];
        }
        b[n] = b[n]/ad[n];
        for(int k=n-1;k>=1;k--){
            b[k] -= au[k]*b[k+1];
            if(k+1 < n && k+1 > n) 
                b[k] = b[k] - ar[k] * b[n];
            b[k] /= ad[k];
        }
    }
}