/** ２次元波動方程式（円領域・陽解法）
 *
 *  <applet code = Wave2d 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 extends Applet{
    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 = 20;
        double tau = 0.005;
        double Tmax = 1.0;
        double dt = 0.05;
        double skip = dt/tau+0.5;
        double hr = 1.0/nr;
        double ht = 2.0*Math.PI/nt;
        double lambdar = tau / hr;
        double lambdat = tau / ht;
        double lambdar2 = lambdar * lambdar;
        double lambdat2 = lambdat * lambdat;
        double Kmax = (Tmax+0.1*tau)/tau;
 
        double[][] u1 = new double[nr+1][nt+1];
        double[][] u2 = new double[nr+1][nt+1];
        double[][] u3 = new double[nr+1][nt+1];

        //初期値代入
        for(int i=0;i<nr;i++){
            for(int j=0;j<nt;j++){
                u1[i][j] = phi(i*hr,j*ht);
            }
        }
        {
            double mean = 0.0;
            for(int j=0;j<nt;j++)
                mean += u1[1][j];
            mean /= nt;
            for(int j=0;j<nt;j++)
                u1[0][j]=mean;
        }
        for(int i=0;i<=nr;i++){
            u1[i][nt]=u1[i][0];
        }
        for(int j=0;j<nt;j++){
            u1[nr][j]=0.0;
        }

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

        for(int i=1;i<nr;i++){
            double ri=i*hr;
            for(int j=0;j<nt;j++){
                int jp1=j+1;
                if(jp1==nt){jp1=0;}
                int jm1=j-1;
                if(jm1==-1){jm1=nt-1;}
                u2[i][j]=(1.0-lambdar2-(lambdat2/(ri*ri)))*u1[i][j]
                    +0.5*(lambdar2*(u1[i-1][j]+u1[i+1][j])
                          +((tau*tau)/(2.0*ri*hr))*(u1[i+1][j]-u1[i-1][j])
                          +(lambdat2/(ri*ri))*(u1[i][jm1]+u1[i][jp1]))
                    +tau*psi(i*hr,j*ht);
            }
        }
        {
            double mean = 0.0;
            for(int j=0;j<nt;j++)
                mean += u2[1][j];
            mean /= nt;
            for(int j=0;j<nt;j++)
                u2[0][j]=mean;
            
        }
        for(int i=0;i<=nr;i++){
            u2[i][nt]=u2[i][0];
        }
        for(int j=0;j<nt;j++){
            u2[nr][j]=0.0;
        }

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

        for(int k=1;k<=Kmax;k++){
            for(int i=1;i<nr;i++){
                double ri=i*hr;
                for(int j=0;j<nt;j++){
                    int jp1=j+1;
                    if(jp1==nt){jp1=0;}
                    int jm1=j-1;
                    if(jm1==-1){jm1=nt-1;}
                    u3[i][j]=2.0*(1.0-lambdar2-(lambdat2/(ri*ri)))*u2[i][j]
                        +lambdar2*(u2[i+1][j]+u2[i-1][j])
                        +((tau*tau)/(2.0*ri*hr))*(u2[i+1][j]-u2[i-1][j])
                        +(lambdat2/(ri*ri))*(u2[i][jp1]+u2[i][jm1])-u1[i][j];
                }
            }

            //r=0
            {
                double mean=0.0;
                for(int j=0;j<nt;j++)
                    mean += u3[1][j];
                mean /= nt;
                for(int j=0;j<nt;j++)
                    u3[0][j]=mean;
            }

            //Dir(r=1)
            for(int j=0;j<nt;j++){
                u3[nr][j]=0.0;
            }
            //θ=2πのところはθ=0のところのコピー
            for(int i=0;i<=nr;i++){
                u3[i][nt]=u3[i][0];
            }

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

            //u1にu2の値をu2にu3の値を代入
            for(int i=0;i<=nr;i++){
                for(int j=0;j<=nt;j++){
                    u1[i][j] = u2[i][j];
                    u2[i][j] = u3[i][j];
                }
            }

            double norm = maxnorm(u3,nr,nt);
            if(norm > 10.0)
                break;

            double t = tau*k;

            //誤差を測る
            double maxe = 0.0;
            double x = 0.0;
            double y = 0.0;
            double e = 0.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) - u3[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 = k*tau;
            }
        }
    }



    //初期条件　ノルム
    public double maxnorm(double[][] u,int nr,int nt){
        int i0 = 0;
        int j0 = 0;
        double absu;
        double tmpmax = Math.abs(u[0][0]);
        for(int i=0;i<=nr;i++){
            for(int j=0;j<=nt;j++){
                if((absu = Math.abs(u[i][j]))>tmpmax){
                    tmpmax = absu;
                    i0 = i;
                    j0 = j;
                }
            }
        }
        return tmpmax;
    }

    double mu01 = 2.404825557695773;
    double mu11 = 3.831705970207512;
    double mu02 = 5.52008;

    //初期条件　φ
    public double phi(double r, double theta)
    {   
        return Bessel.j0(mu01*r)+Bessel.j0(mu02*r);
    }

    //初期条件ψ
    public double psi(double r,double theta)
    {
        return 0.0;
    }

    //厳密解
    public double exactu(double r,double theta,double t){
        return 1.0*Bessel.j0(mu01*r)*Math.cos(mu01*t)
            +2.0*Bessel.j0(mu02*r)*Math.cos(mu02*t);
    } 
}