/**　１次元波動方程式と２次元波動方程式を求めるプログラミングです。
 *   マルチスレッドプログラミングを使っています。
 *   ダブルバッファリングを使っています。
 *   
 *   これをコンパイルするにはMitsuiWorldが必要です。
 *   MitsuiWorldが使える状態であれば普通のコンパイル javac WaveAll_2_0.java
 *   と普通の起動方法 appletviewer Wave1d.java 
 *   でＯＫ。
 *   ただし appletviewer WaveAll_2_0.java としたければコメント欄のどこかに
 *   
 *   <applet code = WaveAll_2_0 width=500 height=500></applet>
 *   
 *   というものが入ってないといけない。
 *   それとMitsuiWorldの使い方によっては
 *
 *   import Mitsui.*;
 *
 *   を消すこと。
 *   MitsuiWorldをパッケージとして使わない人はこれを消してください
 *
 */


import java.applet.*;
import java.awt.*;
import Mitsui.*;
import java.awt.event.*;
import javax.swing.*;
import javax.swing.border.*;
import java.util.*;

public class WaveAll_2_0 extends Applet 
    implements Runnable,ActionListener,ItemListener{

    Thread th = null;                       //スレッド
    Button b_start,b_stop,b_next,b_back,b_reverse;    //各ボタン
    Button b_1dset,b_2dset;
    Label l_time,l_area;                              //各ラベル
    Label2 l_func;
    Choice c_dimension,c_boundary;          //関数を選ぶGUI
    boolean runmove=false;                  //スレッドが動いてるか止まってるか
    boolean reverseflag=false;              //reverseの時に時間をマイナスにする
    
    int im_t,re_t;                          //時間の経過
    int wmax,hmax;                          //画面のサイズ

    //もう一つの画面の定義
    Graphics bg;
    Image buf;

    //画面作成に必要なもの
    MitsuiWorld_2 m; 
    //２次元
    double xmin2,xmax2,ymin2,ymax2,zmin2,zmax2;
    int nx2,ny2;
    double dx2,dy2;
    double tau2;
    //１次元
    double xmin1,xmax1,ymin1,ymax1;
    int nx1;
    double dx1;
    double tau1; 

    int n1=100;
    int n2=50;
    
    double[][] u2_1 = new double[n2+1][n2+1];
    double[][] u2_2 = new double[n2+1][n2+1];
    double[][] u2_3 = new double[n2+1][n2+1];
    double[]   u1_1 = new double[n1+1];
    double[]   u1_2 = new double[n1+1];
    double[]   u1_3 = new double[n1+1];
   
    int nfunc1=0;                            //１次元の関数を変える
    int nfunc2=0;                            //２次元の関数を変える
    int dimension=0;                         //次元を変える/0=１次元,1=２次元
    int boundary=0;                          //境界条件を変える

    //視覚を変える変数
    int lx,ly,lz;                            //各軸の長さ
    int angx,angy;                           //視覚の角度
    int x0,y0;                               //原点の位置

    //初期設定    
    public void init(){
        setLayout(new BorderLayout());
        Panel p  = new Panel();            //NORTH地区を枠決めしてる
        p.setLayout(new GridLayout(3,1,0,0));
        Panel p1 = new Panel();            //１段目
        Panel p2 = new Panel();            //２段目
        Panel p3 = new Panel();            //３段目
        

        //startボタンの設置
        b_start=new Button("START");
        b_start.addActionListener(this);
        p2.add(b_start);
        
        //stopボタンの設置
        b_stop=new Button("STOP");
        b_stop.addActionListener(this);
        p2.add(b_stop);

        //reverseボタンの設置
        b_reverse=new Button("REVERSE");
        b_reverse.addActionListener(this);
        p2.add(b_reverse);

        //backボタンの設置
        b_back=new Button("BACK");
        b_back.addActionListener(this);
        p2.add(b_back);

        //nextボタンの設置
        b_next=new Button("NEXT");
        b_next.addActionListener(this);
        p2.add(b_next);

        //dimensionチョイスの設置
        c_dimension=new Choice();
        c_dimension.addItem("１次元波動方程式");
        c_dimension.addItem("２次元波動方程式");
        c_dimension.addItemListener(this);
        p1.add(c_dimension);

        //boundaryチョイスの設置
        c_boundary=new Choice();
        c_boundary.addItem("Dirichlet");
        c_boundary.addItem("Neumann");
        c_boundary.addItem("仮想格子点");
        c_boundary.addItemListener(this);
        p1.add(c_boundary);

        //１次元のセッティングを行うボタン
        b_1dset = new Button("1Dsetting");
        b_1dset.addActionListener(this);
        p1.add(b_1dset);

        //２次元のセッティングを行うボタン
        b_2dset = new Button("2Dsetting");
        b_2dset.addActionListener(this);
        p1.add(b_2dset);
        
        //ラベルl_timeの設置　時間を表示する。
        l_time = new Label("           ");
        p2.add(l_time);

        //ラベルl_areaの設置　描写範囲の表示
        l_area = new Label("                                                 ");
        p3.add(l_area);

        //ラベルl_funcの設置　選んだ関数の表示
        l_func = new Label2("φ=sin(πx)　　　　　　　　　　　　   　　",
                            "ψ=0　　　　　　　　　　　　　    　　　　");
        p3.add(l_func);

        p.add(p1);
        p.add(p2);
        p.add(p3);
        add(p,"North");
        
        //画面の横と縦の長さを得る。
        wmax = getSize().width;
        hmax = getSize().height;

        //描写範囲を決める
        //２次元[xmin2,xmax2]×[ymin2,ymax2]×[zmin2,zmax2]
        xmin2 = -2.0;
        xmax2 = 2.5;
        ymin2 = -2.0;
        ymax2 = 2.5;
        zmin2 = -1.0;
        zmax2 = 1.0;

        //１次元[xmin1,xmax1]×[ymin1,ymax1]
        xmin1 = -0.2;
        xmax1 = 1.2;
        ymin1 = -1.2;
        ymax1 = 1.8;

        //作図用の描写画面の作成
        buf = createImage(wmax,hmax);
        bg = buf.getGraphics();
        
        //MitsuiWorld
        m = new MitsuiWorld_2();
        m.setGraphics(bg);
        m.setScreenSize(wmax,hmax);
                        
        //分割数
        nx1 = 100;                  //1次元のｘ分割数
        nx2 = 40;                   //２次元のxの分割数大きすぎると動作が遅い
        ny2 = 40;                   //２次元のyの分割数大きすぎると動作が遅い
        
        //時刻間隔
        tau1 = 0.01;
        tau2 = 0.07;

        //視覚を変える変数の初期値
        lx=0;
        ly=0;
        lz=0;
        x0=y0=0;
        angx=angy=0;
    }
    
    //起動と同時にスレッドを走らせる
    public void start(){
        if(th==null){
            th=new Thread(this);
            th.start();
        }
    }
    
    //スレッドの実装
    public void run(){
        while(true){                 //起動中スレッドを走らせる。
            while(!runmove){           //startボタンを押すことによって
            }                           //ここからぬけ、スレッドが有効になる。
            if(!reverseflag)
                re_t++;
            else
                re_t--;
            
            im_t++;

            bg.clearRect(0,0,wmax,hmax);        //画面をクリア
            if(dimension==0){                   //１次元波動方程式のとき
                l_time.setText("t="+re_t*tau1);
                double[] u = f1(im_t);
                m.setColor(Color.black);
                drawAxis();
                m.setColor(Color.pink);
                m.move(0.0,u[0]);
                for(int i=1;i<=nx1;i++)
                    m.draw(i*dx1,u[i]);
            }
            if(dimension==1){                   //２次元波動方程式のとき
                l_time.setText("t="+re_t*tau2);
                double[][] u = f2(im_t);           //計算した結果を格納する
                m.hideBirdView(u,nx2,ny2,1);         //メモリ内で描く
            }
            repaint();
            try{
                Thread.sleep(100);          //100ミリ秒ストップ
            }catch(InterruptedException e){}
        }
    }
    
    //終了と同時にスレッドも終了
    public void stop(){
        if(th!=null){
            th = null;
        }
    }
    
    public void actionPerformed(ActionEvent e){
        //startボタンを押したときのアクション
        if(e.getSource()==b_start){
            runmove=true;                //再びグラフを描き始める
        }

        //stopボタンを押したときのアクション
        if(e.getSource()==b_stop){
            if(!runmove){                          //一時停止のとき
                reverseflag=false;                 //初期化
                im_t=0;
                re_t=0;
                repaint();
            }
            else{                                 //これで一時停止 
                runmove=false;
            }
        }

        //backボタンを押したときのアクション
        if(e.getSource()==b_back){
            if(!reverseflag)
                re_t--;
            else
                re_t++;

            if(!runmove){
                if(dimension==0){
                    for(int i=0;i<=nx1;i++){
                        u1_2[i] = u1_1[i];
                        u1_1[i] = u1_3[i];
                    }
                    double[] u = f1(im_t);
                    for(int i=0;i<=nx1;i++){
                        u1_2[i] = u1_1[i];
                        u1_1[i] = u1_3[i];
                        u1_3[i] = u1_2[i];
                    }

                    l_time.setText("t="+(double)re_t*tau1);
                    bg.clearRect(0,0,wmax,hmax);        //画面をクリア 
                    m.setColor(Color.black);
                    drawAxis();
                    m.setColor(Color.pink);
                    m.move(0.0,u1_3[0]);
                    for(int i=1;i<=nx1;i++)
                        m.draw(i*dx1,u1_3[i]);
                }
                if(dimension==1){
                    for(int i=0;i<=nx2;i++){                
                        for(int j=0;j<=ny2;j++){
                            u2_2[i][j] = u2_1[i][j];
                            u2_1[i][j] = u2_3[i][j];
                        }
                    }
                    double[][] u = f2(im_t);
                    for(int i=0;i<=nx2;i++){                
                        for(int j=0;j<=ny2;j++){
                            u2_2[i][j] = u2_1[i][j];
                            u2_1[i][j] = u2_3[i][j];
                            u2_3[i][j] = u2_2[i][j];
                        }
                    }
                    
                    l_time.setText("t="+(double)re_t*tau2);
                    bg.clearRect(0,0,wmax,hmax);        //画面をクリア
                    m.hideBirdView(u2_3,nx2,ny2,1);         //メモリ内で描く
                }
                repaint();              
            }
        }

        //nextボタンを押したときのアクション
        if(e.getSource()==b_next){
            if(!runmove){                      //一時停止状態のときのみ
                if(!reverseflag)
                re_t++;
            else
                re_t--;                          //次の時間のグラフを見れる。

                bg.clearRect(0,0,wmax,hmax);        //画面をクリア
                if(dimension==0){
                    l_time.setText("t="+(double)re_t*tau1);
                    double[] u = f1(im_t);
                    m.setColor(Color.black);
                    drawAxis();
                    m.setColor(Color.pink);
                    m.move(0.0,u[0]);
                    for(int i=1;i<=nx1;i++)
                        m.draw(i*dx1,u[i]);
                }
                if(dimension==1){
                    l_time.setText("t="+(double)re_t*tau2);
                    double[][] u = f2(im_t);           //計算した結果を格納する
                    m.hideBirdView(u,nx2,ny2,1);         //メモリ内で描く
                }
                repaint();              
            }
        }

        //reverseボタンを押したときのアクション
        if(e.getSource()==b_reverse){
            //逆流してるかしてないか。
            if(reverseflag)
                reverseflag=false;
            if(!reverseflag)
                reverseflag=true;

            if(dimension==0){
                double[]u = f1(im_t);
                for(int i=0;i<=nx1;i++){
                    u1_2[i] = u1_1[i];
                    u1_1[i] = u1_3[i];
                    u1_3[i] = u1_2[i];
                }
            }
            if(dimension==1){
                double[][] u = f2(im_t); 
                for(int i=0;i<=nx2;i++){
                    for(int j=0;j<=ny2;j++){
                        u2_2[i][j]=u2_1[i][j];
                        u2_1[i][j]=u2_3[i][j];
                        u2_3[i][j]=u2_2[i][j];
                    }
                }
            }
        }

        //１次元波動方程式の設定をするダイアログを作っている。
        if(e.getSource()==b_1dset){
            String s = ("初期条件を選んでください");

            //初期値選択欄の作成
            Choice c_func1 = new Choice();
            c_func1.addItem("ひたすら振動を繰り返す波");
            c_func1.addItem("別れては重なる波");
            c_func1.addItem("右方向に流れていく波");
            c_func1.addItem("二つの波の合体");
            
            //分割数記入欄の作成
            NumericField inputN = new NumericField(""+nx1);
            inputN.setBorder
                (new TitledBorder("分割数を指定してください(=<100)"));

            //時間刻み記入欄の作成
            NumericField  inputTau = new NumericField(""+tau1);
            inputTau.setBorder
                (new TitledBorder("時間刻みを指定してください"));

            //描写範囲指定欄の作成
            NumericField  inputArea =
                new NumericField(xmin1+" "+xmax1+"         "+ymin1+" "+ymax1);
            inputArea.setBorder
                (new TitledBorder("描写範囲を指定してください　数の間には空白を入れて下さい。"));

            //オブジェクトを一つにまとめて
            Object[] obj = {s,c_func1,inputN,inputTau,inputArea};
           
            //ダイアログの作成
            int ans = JOptionPane.showConfirmDialog(this,obj,
                                              "１次元波動方程式の設定",
                                               JOptionPane.OK_CANCEL_OPTION);

            if(ans==0){                           //OKを押した時
                nfunc1  = c_func1.getSelectedIndex();

                if(dimension==0){
                    //式の記述
                    switch(nfunc1){
                    case 0:
                        l_func.setText("φ=sin(πx)","ψ=0");
                        break;
                    case 1:
                        l_func.setText("φ=8x-3(3/8<x<=1/2),-8x+5(1/2<x<=5/8)",
                                       "ψ=0");
                        break;
                    case 2:
                        l_func.setText("φ=(sin(π(10x-0.5))+1)/4",
                                       "ψ=-2.5πcos(π(10x-0.5))");
                        break;
                    case 3:
                        l_func.setText("φ1=(sin(π(10x-0.5))+1)/4",
                                       "φ2=(sin(π(20(x-0.9)-0.5))+1)/8");
                        break;
                    default:
                        l_func.setText("未定","未定");
                        break;
                    }
                }
                
                nx1 = Integer.parseInt(inputN.getText().trim());
                tau1 = Double.valueOf(inputTau.getText().trim()).doubleValue();
                
                //描写範囲を読み込む時は「分割」をつかう。
                String str = inputArea.getText().trim();
                //空白区切りで分割する
                StringTokenizer st = new StringTokenizer(str," ");
                if(st.countTokens()==4){             //分割数が４のとき
                    xmin1 = Double.valueOf(st.nextToken()).doubleValue();
                    xmax1 = Double.valueOf(st.nextToken()).doubleValue();
                    ymin1 = Double.valueOf(st.nextToken()).doubleValue();
                    ymax1 = Double.valueOf(st.nextToken()).doubleValue();
                }

                                
                runmove  = false;                          //初期状態に戻す
                reverseflag=false;
                im_t=0;
                re_t=0;
                repaint();
            }
        }

        if(e.getSource()==b_2dset){
            String s = ("初期条件を選んでください");

            Choice c_func2 = new Choice();
            c_func2.addItem("φ=sin(πx)+sin(πy),ψ=0");
            c_func2.addItem("定常波");
            c_func2.addItem("平面波");
            c_func2.addItem("平面波の合体");
            c_func2.addItem("変な波");
            c_func2.addItem("変な波右方向");
            
            //ダイアログは持ち込むオブジェクトを縦に並べる。しかし
            //パネルを使うことにより横にもコンテンツを増やした。
            Panel pdivide = new Panel();        
            JLabel N = new JLabel("分割数");
            NumericField inputNx = new NumericField(""+nx2,10);
            NumericField inputNy = new NumericField(""+ny2,10);
            inputNx.setBorder(new TitledBorder("Nx(=<50)"));
            inputNy.setBorder(new TitledBorder("Ny(=<50)"));
            pdivide.add(N);
            pdivide.add(inputNx);
            pdivide.add(inputNy);
           
            NumericField  inputTau = new NumericField(""+tau2);
            inputTau.setBorder
                (new TitledBorder("時間刻みを指定してください"));

            NumericField  inputArea =
                new NumericField(xmin2+" "+xmax2+"         "+
                                 ymin2+" "+ymax2+"         "+zmin2+" "+zmax2);
            inputArea.setBorder
                (new TitledBorder("描写範囲の指定　数の間には空白を入れて下さい。"));
            //角度を変える欄の作成
            Panel pangle = new Panel();        
            JLabel A = new JLabel("角度変え（ラジアン）");
            NumericField inputAngx = new NumericField(""+angx,10);
            NumericField inputAngy = new NumericField(""+angy,10);
            inputAngx.setBorder(new TitledBorder("水平方向"));
            inputAngy.setBorder(new TitledBorder("垂直方向"));
            pangle.add(A);
            pangle.add(inputAngx);
            pangle.add(inputAngy);
            
            //原点移動の欄の作成
            Panel porigin = new Panel();        
            JLabel O = new JLabel("原点移動（ピクセル）");
            NumericField inputX0 = new NumericField(""+x0,10);
            NumericField inputY0 = new NumericField(""+y0,10);
            inputX0.setBorder(new TitledBorder("右向き"));
            inputY0.setBorder(new TitledBorder("上向き"));
            porigin.add(O);
            porigin.add(inputX0);
            porigin.add(inputY0);

            //軸の長さ変更欄の作成
            Panel plength = new Panel();        
            JLabel L = new JLabel("サイズ変え");
            NumericField inputLx = new NumericField(""+lx,8);
            NumericField inputLy = new NumericField(""+ly,8);
            NumericField inputLz = new NumericField(""+lz,8);
            inputLx.setBorder(new TitledBorder("ｘ軸"));
            inputLy.setBorder(new TitledBorder("ｙ軸"));
            inputLz.setBorder(new TitledBorder("ｚ軸"));
            plength.add(L);
            plength.add(inputLx);
            plength.add(inputLy);
            plength.add(inputLz);

            Object[] obj = {s,c_func2,pdivide,inputTau,inputArea,pangle,
                            porigin,plength};
           
            int ans = JOptionPane.showConfirmDialog(this,obj,
                                              "２次元波動方程式の設定",
                                               JOptionPane.OK_CANCEL_OPTION);
            if(ans==0){
                nfunc2  = c_func2.getSelectedIndex();

                if(dimension==1){
                    //式の記述
                    switch(nfunc2){
                    case 0:
                        l_func.setText("φ=sin(πx)+sin(πy)",
                                       "ψ=0");
                        break;
                    case 1:
                        l_func.setText("m=(A+B)/2,n=(C+D)/2,r=(B-A)/5",
                                       "φ=sqrt(r^2-(x-m)^2-(x-n)^2)(m^2+n^2=<r^2)");
                        break;
                    case 2:
                        l_func.setText("φ=(sin(π(2x-0.5))+1)/4",
                                       "ψ=-2πcos(π(2x-0.5))/4");
                        break;
                    case 3:
                        l_func.setText("φ1=(sin(π(2x-0.5))+1)/4",
                                       "φ2=(sin(π(2(x-B+1)-0.5))+1)/6");
                        break;
                    default:
                        l_func.setText("未定","未定");
                        break;
                    }
                }

                nx2 = Integer.parseInt(inputNx.getText().trim());
                ny2 = Integer.parseInt(inputNy.getText().trim());
                
                tau2 = Double.valueOf(inputTau.getText().trim()).doubleValue();

                String str = inputArea.getText().trim();
                StringTokenizer st = new StringTokenizer(str," ");
                if(st.countTokens()==6){
                    xmin2 = Double.valueOf(st.nextToken()).doubleValue();
                    xmax2 = Double.valueOf(st.nextToken()).doubleValue();
                    ymin2 = Double.valueOf(st.nextToken()).doubleValue();
                    ymax2 = Double.valueOf(st.nextToken()).doubleValue();
                    zmin2 = Double.valueOf(st.nextToken()).doubleValue();
                    zmax2 = Double.valueOf(st.nextToken()).doubleValue();
                }

                angx = Integer.parseInt(inputAngx.getText().trim());
                angy = Integer.parseInt(inputAngy.getText().trim());

                x0 = Integer.parseInt(inputX0.getText().trim());
                y0 = Integer.parseInt(inputY0.getText().trim());

                lx = Integer.parseInt(inputLx.getText().trim());
                ly = Integer.parseInt(inputLy.getText().trim());
                lz = Integer.parseInt(inputLz.getText().trim());

                runmove  = false;                          //初期状態に戻す
                reverseflag=false;
                im_t=0;
                re_t=0;
                repaint();
            }
        }
    }

    //dimension,boundary,fundチョイスを選んだとき
    public void itemStateChanged(ItemEvent e){
        dimension=c_dimension.getSelectedIndex();  //次元番号を変える
        boundary =c_boundary.getSelectedIndex();   //境界条件番号を変える
        runmove  = false;                          //初期状態に戻す
        reverseflag=false;
        im_t=0;
        re_t=0;
        repaint();
    }
        
    public void paint(Graphics g){
        if(im_t==0){ 
            bg.clearRect(0,0,wmax,hmax);
            if(dimension==0){
                dx1  = 1.0/nx1;             //１次元のｘ格子点
                l_time.setText("t=0.0");    //時刻と描写範囲の表示
                l_area.setText("["+xmin1+","+xmax1+"]×["+ymin1+","+ymax1+"]");
                m.setArea(xmin1,xmax1,ymin1,ymax1);

                double[] u = f1(im_t);
                m.setColor(Color.black);
                drawAxis();
                m.setColor(Color.pink);
                m.move(0.0,u[0]);
                for(int i=1;i<=nx1;i++)
                    m.draw(i*dx1,u[i]);
            }
            if(dimension==1){
                dx2 = (xmax2-xmin2)/nx2;    //２次元のｘ格子点
                dy2 = (ymax2-ymin2)/ny2;    //２次元のｙ格子点
                l_time.setText("t=0.0");    //時刻と描写範囲の表示
                l_area.setText("["+xmin2+","+xmax2+"]×["+ymin2+","+ymax2+
                               "]×["+zmin2+","+zmax2+"]");
                //視覚を変えている
                m.changeAngle(angx,angy);
                m.changePosition(x0,y0);
                m.changeSize(lx,ly,lz);

                //描き始める。
                m.setArea2(xmin2,xmax2,ymin2,ymax2,zmin2,zmax2);
                double[][] u = f2(im_t);
                m.setColor(Color.pink);
                m.hideBirdView(u,nx2,ny2,1);
            }
        }
        g.drawImage(buf,0,0,this);            //もう一つの画面をくっつける
    }

    public double[][] f2(int t){
        double lambdax = tau2 / dx2;
        double lambday = tau2 / dy2;
        double lambdax2 = lambdax * lambdax;
        double lambday2 = lambday * lambday;
        
        if(t==0){                                 //t=0のとき
            for(int i=0;i<=nx2;i++)
                for(int j=0;j<=ny2;j++)
                    u2_1[i][j] = phi(xmin2+i*dx2,ymin2+j*dy2);     //初期値代入
            
            return u2_1;
        }
        
        if(t==1){  
            for(int i=1;i<nx2;i++)
                for(int j=1;j<ny2;j++)
                    u2_2[i][j] = (1.0-lambdax2-lambday2) * u2_1[i][j] + 
                        0.5 * lambdax2 * (u2_1[i-1][j]+u2_1[i+1][j]) +
                        0.5 * lambday2 * (u2_1[i][j-1]+u2_1[i][j+1]) +
                        tau2 * psi(xmin2+i*dx2,ymin2+j*dy2);

            if(boundary==0){                          //Dirichlet境界条件
                for(int i=0;i<=nx2;i++)
                    u2_2[i][0]=u2_2[i][ny2]=0.0;
                for(int i=0;i<=ny2;i++)
                    u2_2[0][i]=u2_2[nx2][i]=0.0;
            }
            if(boundary==1){                          //Neumann境界条件
                for(int i=0;i<=nx2;i++){
                    u2_2[i][0]=u2_2[i][1];
                    u2_2[i][ny2]=u2_2[i][ny2-1];
                }
                for(int i=0;i<=ny2;i++){
                    u2_2[0][i]=u2_2[1][i];
                    u2_2[nx2][i]=u2_2[nx2-1][i];
                }
            }
            if(boundary==2){                         //仮想格子点
                for(int i=1;i<nx2;i++){
                    u2_2[i][0]=(1.0-lambdax2-lambday2) * u2_1[i][0] + 
                        0.5 * lambdax2 * (u2_1[i-1][0]+u2_1[i+1][0]) +
                        0.5 * lambday2 * (2*u2_1[i][1]) +
                        tau2 * psi(xmin2+i*dx2,ymin2);
                    u2_2[i][ny2]=(1.0-lambdax2-lambday2) * u2_1[i][ny2] + 
                        0.5 * lambdax2 * (u2_1[i-1][ny2]+u2_1[i+1][ny2]) +
                        0.5 * lambday2 * (2*u2_1[i][ny2-1]) +
                        tau2 * psi(xmin2+i*dx2,ymin2+ny2*dy2);
                }
                for(int j=1;j<ny2;j++){
                    u2_2[0][j]=(1.0-lambdax2-lambday2) * u2_1[0][j] + 
                        0.5 * lambdax2 * (2*u2_1[1][j]) +
                        0.5 * lambday2 * (u2_1[0][j-1]+u2_1[0][j+1]) +
                        tau2 * psi(xmin2,ymin2+j*dy2);
                    u2_2[nx2][j]=(1.0-lambdax2-lambday2) * u2_1[nx2][j] + 
                        0.5 * lambdax2 * (2*u2_1[nx2-1][j]) +
                        0.5 * lambday2 * (u2_1[nx2][j-1]+u2_1[nx2][j+1]) +
                        tau2 * psi(xmin2+nx2*dx2,ymin2+j*dy2);
                }
                u2_2[0][0]=(1.0-lambdax2-lambday2) * u2_1[0][0] + 
                        0.5 * lambdax2 * (2*u2_1[1][0]) +
                        0.5 * lambday2 * (2*u2_1[0][1]) +
                        tau2 * psi(xmin2,ymin2);
                u2_2[nx2][0]=(1.0-lambdax2-lambday2) * u2_1[nx2][0] + 
                        0.5 * lambdax2 * (2*u2_1[nx2-1][0]) +
                        0.5 * lambday2 * (2*u2_1[nx2][1]) +
                        tau2 * psi(xmin2+nx2*dx2,ymin2);
                u2_2[0][ny2]=(1.0-lambdax2-lambday2) * u2_1[0][ny2] + 
                        0.5 * lambdax2 * (2*u2_1[1][ny2]) +
                        0.5 * lambday2 * (u2_1[0][ny2-1]) +
                        tau2 * psi(xmin2,ymin2+ny2*dy2);
                u2_2[nx2][ny2]=(1.0-lambdax2-lambday2) * u2_1[nx2][ny2] + 
                        0.5 * lambdax2 * (2*u2_1[nx2-1][ny2]) +
                        0.5 * lambday2 * (2*u2_1[nx2][ny2-1]) +
                        tau2 * psi(xmin2+nx2*dx2,ymin2+ny2*dy2);
                
            }
            return u2_2;
        }
        
        else{                                     //t>=2のとき
            for(int i=1;i<nx2;i++)
                for(int j=1;j<ny2;j++)
                    u2_3[i][j] = 2 * (1.0-lambdax2 -lambday2) * u2_2[i][j] + 
                        lambdax2 * (u2_2[i-1][j]+u2_2[i+1][j]) +
                        lambday2 * (u2_2[i][j-1]+u2_2[i][j+1]) - u2_1[i][j];
            
            if(boundary==0){                          //Dirichlet境界条件
                for(int i=0;i<=nx2;i++)
                    u2_3[i][0]=u2_3[i][ny2]=0.0;
                for(int i=0;i<=ny2;i++)
                    u2_3[0][i]=u2_3[nx2][i]=0.0;
            }
            if(boundary==1){                          //Neumann境界条件
                for(int i=0;i<=nx2;i++){
                    u2_3[i][0] = u2_3[i][1] ;
                    u2_3[i][ny2]= u2_3[i][ny2-1];
                }
                for(int i=0;i<=ny2;i++){
                    u2_3[0][i] = u2_3[1][i];
                    u2_3[nx2][i]= u2_3[nx2-1][i];
                }
            }
            if(boundary==2){                         //仮想格子点
                for(int i=1;i<nx2;i++){
                    u2_3[i][0] = 2 * (1.0-lambdax2 -lambday2) * u2_2[i][0] + 
                        lambdax2 * (u2_2[i-1][0]+u2_2[i+1][0]) +
                        lambday2 * (2*u2_2[i][1]) - u2_1[i][0];
                    u2_3[i][ny2]=2 * (1.0-lambdax2 -lambday2) * u2_2[i][ny2] + 
                        lambdax2 * (u2_2[i-1][ny2]+u2_2[i+1][ny2]) +
                        lambday2 * (2*u2_2[i][ny2-1]) - u2_1[i][ny2];
                }
                for(int j=1;j<ny2;j++){
                    u2_3[0][j] = 2 * (1.0-lambdax2 -lambday2) * u2_2[0][j] + 
                        lambdax2 * (2*u2_2[1][j]) +
                        lambday2 * (u2_2[0][j-1]+u2_2[0][j+1]) - u2_1[0][j];
                    u2_3[nx2][j]=2 * (1.0-lambdax2 -lambday2) * u2_2[nx2][j] + 
                        lambdax2 * (2*u2_2[nx2-1][j]) +
                        lambday2 *(u2_2[nx2][j-1]+u2_2[nx2][j+1])-u2_1[nx2][j];
                }
                u2_3[0][0] = 2 * (1.0-lambdax2 -lambday2) * u2_2[0][0] + 
                        lambdax2 * (2*u2_2[1][0]) +
                        lambday2 * (2*u2_2[0][1]) - u2_1[0][0];
                u2_3[nx2][0]= 2 * (1.0-lambdax2 -lambday2) * u2_2[nx2][0] + 
                        lambdax2 * (2*u2_2[nx2-1][0]) +
                        lambday2 * (2*u2_2[nx2][1]) - u2_1[nx2][0];
                u2_3[0][ny2]= 2 * (1.0-lambdax2 -lambday2) * u2_2[0][ny2] + 
                        lambdax2 * (2*u2_2[1][ny2]) +
                        lambday2 * (2*u2_2[0][ny2-1]) - u2_1[0][ny2];
                u2_3[nx2][ny2]=2 * (1.0-lambdax2 -lambday2) * u2_2[nx2][ny2] + 
                        lambdax2 * (2*u2_2[nx2-1][ny2]) +
                        lambday2 * (2*u2_2[nx2][ny2-1]) - u2_1[nx2][ny2];
            }
            
            //u1にu2の値をu2にu3の値を代入
            for(int i=0;i<=nx2;i++){
                for(int j=0;j<=ny2;j++){
                    u2_1[i][j] = u2_2[i][j];
                    u2_2[i][j] = u2_3[i][j];
                }
            }
            
            return u2_3;
        }
        
    }
    
    public double[] f1(int t){
        double lambda = tau1/dx1;
        double lambda2 = lambda * lambda;
        
        if(t==0){
            for(int i=0;i<=nx1;i++)
                u1_1[i] = phi(i*dx1); 
            return u1_1;
        }
        
        if(t==1){
            for(int i=1;i<nx1;i++)
                u1_2[i] = (1-lambda2) * u1_1[i] + 
                    0.5 * lambda2 * (u1_1[i-1]+u1_1[i+1]) + tau1 * psi(i*dx1);
            
            if(boundary==0){                       //Dirichlet境界条件
                u1_2[0]=u1_2[nx1]=0;
            }
            if(boundary==1){                       //Neumann境界条件
                u1_2[0]=u1_2[1];
                u1_2[nx1]=u1_2[nx1-1];
            }
            if(boundary==2){                       //仮想格子点
                u1_2[0] = (1-lambda2) * u1_1[0] + 
                    0.5 * lambda2 * (2*u1_1[1]) + tau1 * psi(0.0);
                u1_2[nx1] = (1-lambda2) * u1_1[nx1] + 
                    0.5 * lambda2 * (2*u1_1[nx1-1]) + tau1 * psi(nx1*dx1);
            }
            
            return u1_2;
        }
        
        else{
            for(int i=1;i<nx1;i++)
                u1_3[i] = 2 * (1.0-lambda2) * u1_2[i] + 
                    lambda2 * (u1_2[i-1]+u1_2[i+1]) - u1_1[i];
            
            if(boundary==0){                       //Dirichlet境界条件
                u1_3[0] = u1_3[nx1] = 0;
            }
            if(boundary==1){                       //Neumann境界条件 
                u1_3[0]=u1_3[1];
                u1_3[nx1]=u1_3[nx1-1];
            }
            if(boundary==2){                       //仮想格子点
                u1_3[0] = 2 * (1.0-lambda2) * u1_2[0] + 
                    lambda2 * (2*u1_2[1]) - u1_1[0];
                u1_3[nx1] = 2 * (1.0-lambda2) * u1_2[nx1] + 
                    lambda2 * (2*u1_2[nx1-1]) - u1_1[nx1];
            }
            
            //u1にu2の値をu2にu3の値を代入
            for(int i=0;i<nx1;i++){
                u1_1[i] = u1_2[i];
                u1_2[i] = u1_3[i];
            }
            
            return u1_3;
        }
    }

    //初期値φ
    public double phi(double x,double y){
        double cx = (xmax2+xmin2)/2;
        double cy = (ymax2+ymin2)/2;
        double r  = (xmax2-xmin2)/5;

        if(nfunc2==0){
            return Math.sin(Math.PI*x)/2+Math.sin(Math.PI*y)/2;
        }
        if(nfunc2==1){
            if((x-cx)*(x-cx)+(y-cy)*(y-cy)<=r*r)
                return Math.sqrt(r*r-(x-cx)*(x-cx)-(y-cy)*(y-cy));
            else
                return 0.0;
        }
        if(nfunc2==2){
            if(xmin2<=x && x<=xmin2+1)
                return (Math.sin(Math.PI*(x*2-0.5))+1)/4;
            else 
                return 0.0;
        }
        if(nfunc2 == 3){
            if(xmin2<=x && x<=xmin2+1)
                return (Math.sin(Math.PI*(x*2-0.5))+1)/4;
            if(xmax2-1<=x && x<=xmax2)
                return (Math.sin(Math.PI*((x-xmax2+1)*2-0.5))+1)/6;
            else
                return 0.0;
        }
        if(nfunc2 == 4){
            if(x*x+y*y<1)
                return ((Math.cos(Math.PI*x)+Math.cos(Math.PI*y))+1)/2;
            else
                return 0.0;
        }
        if(nfunc2 == 5){
            if(x*x+y*y<1)
                return ((Math.cos(Math.PI*x)+Math.cos(Math.PI*y))+1)/2;
            else
                return 0.0;
        }
        else
            return 0.0;
    }
    
    //初期条件ψ
    public double psi(double x,double y){
        double cx = (xmax2+xmin2)/2;
        double cy = (ymax2+ymin2)/2;
        double r  = (xmax2-xmin2)/5;

        if(nfunc2==2){
            if(xmin2<=x && x<=xmin2+1)
                return -2*Math.PI*Math.cos(Math.PI*(x*2-0.5))/4;
            else
                return 0.0;
        }
        if(nfunc2 == 3){
            if(xmin2<=x && x<=xmin2+1)
                return -2*Math.PI*Math.cos(Math.PI*(x*2-0.5))/4;
            if(xmax2-1<=x && x<=xmax2)
                return 2*Math.PI*Math.cos(Math.PI*((x-xmax2+1)*2-0.5))/6;
            else
                return 0.0;
        }
        if(nfunc2 == 4){
            if(x*x+y*y<1)
                return (-Math.PI*Math.sin(Math.PI*x)+
                        -Math.PI*Math.sin(Math.PI*y))/2;
            else
                return 0.0;
        }
        if(nfunc2 ==5){
            if(x*x+y*y<1)
                return Math.PI*Math.sin(Math.PI*x)/2;
            else
                return 0.0;
        }
        else
            return 0.0;
    }
    
    //初期値φ
    public double phi(double x){
        if(nfunc1 == 0){
            return Math.sin(Math.PI*x);
        }
        if(nfunc1 == 1){
            if(0.375<x && x<=0.5)
                return 8.0*x-3.0;
            if(0.5<x && x<=0.625)
                return -8.0*x+5.0;
            else 
                return 0.0;
        }
        if(nfunc1 == 2){
            if(0.0<=x && x<=0.2)
                return (Math.sin(Math.PI*(x*10-0.5))+1)/4;
            else
                return 0.0;
        }
        if(nfunc1 == 3){
            if(0.0<=x && x<=0.2)
                return (Math.sin(Math.PI*(x*10-0.5))+1)/4;
            if(0.9<=x && x<=1.0)
                return (Math.sin(Math.PI*((x-0.9)*20-0.5))+1)/8;
            else
                return 0.0;
        }
        else
            return 0.0;
    }

    //初期条件ψ  
    public double psi(double x){
        if(nfunc1 == 0){
            return 0.0;
        }
        if(nfunc1 == 1){
            return 0.0;
        }
        if(nfunc1 == 2){
            if(0.0<=x && x<=0.2)
                return -2.5*Math.PI*Math.cos(Math.PI*(10*x-0.5));
            else
                return 0.0;
        }
        if(nfunc1 == 3){
            if(0.0<=x && x<=0.2)
                return -2.5*Math.PI*Math.cos(Math.PI*(10*x-0.5));
            if(0.9<=x && x<=1.0)
                return 2.5*Math.PI*Math.cos(Math.PI*((x-0.9)*20-0.5));
            else
                return 0.0;
        }
        else
            return 0.0;
    }
    
    public void drawAxis(){
        m.move(0.0,0.0); m.draw(1.025,0.0);
        for(double i=0.0;i<=1.0;i+=0.1){
            m.move(i,0.025); m.draw(i,-0.025);
        }
        m.move(0.0,1.05); m.draw(0.0,-1.05);
        for(double i=1.0;i>=-1;i-=0.2){
            m.move(-0.01,i); m.draw(0.01,i);
        }
    }

    //数字しか入らないテキストフィールドの作成
    class NumericField extends JTextField{
        public NumericField(String begin){
            super(begin);
            enableEvents(AWTEvent.KEY_EVENT_MASK);
        }

        public NumericField(String begin,int num){
            super(begin,num);
            enableEvents(AWTEvent.KEY_EVENT_MASK);
        }

        protected void processKeyEvent(KeyEvent e){
            String validValues = "0123456789.- ";
            int code = e.getKeyCode();
            switch(code){
            case 39:break;                         //右カーソル
            case 37:break;                         //左カーソル
            default:
                char chr = e.getKeyChar();
                if(chr >= ' ' && validValues.indexOf( chr )==-1){
                    return;
                }
                break;
            }
            super.processKeyEvent(e);
        }
    }

    //２行にまたがるラベル。
    class Label2 extends JComponent{
        JLabel row1,row2;
        
        public Label2(String s1,String s2){
            setLayout(new BoxLayout(this,BoxLayout.Y_AXIS));
            row1 = new JLabel(s1);
            row2 = new JLabel(s2);
            add(row1);
            add(row2);
        }

        public void setText(String s1,String s2){
            row1.setText(s1);
            row2.setText(s2);
        }
    }
}