//<APPLET code="Oregonator_vector_mk.class"width=1050 height=650></APPLET>

import java.applet.*;
import java.awt.*;
import java.awt.event.*;
import corejava.*;

public class Oregonator_vector_mk extends Applet implements MouseListener,ActionListener {

    // 最終時刻 tt の初期値
    static final double init_tt = 20.0;
    private double eps_tol = 1e-5;
    private boolean verbose = false;
    //スクリーン1
    private Image img1;
    private Graphics gra1;
    
    //スクリーン2
    private Image img2;
    private Graphics gra2;

    // RKF45 のための引数
    private double h;
    private double hh;
    private double [] t;
    private double [] x;
    private double [] z;
    private double next_x, next_z, error, xn, zn;

    //ユーザーとのインターフェイス
    private double q,eps,xx0,zz0,tt,tt0,cf, 
        x_max, x_min, y_max, y_min, x_margin, y_margin;
    private int n, maxn;
    private Label la1,la2,la3,la4,la5,la6,la7,la8,la9,la10,la11,la12,la13;
    private TextField tf1,tf2,tf3,tf4,tf5,tf6,tf7,tf8,tf9,tf10,tf11,tf12,tf13;
    private Button bt1,bt2,bt3,bt4,bt5;
    private TextArea ta1, ta2;

    // x-z 座標系の変換パラメーター
    private double ratiox,ratioy,X0,Y0;
    private int WX = 650, WY = 650;

    private boolean IsIn1(double x, double y) {
        return (x_min <= x && x <= x_max && y_min <= y && y <= y_max);
    }

    // x-z 座標変換の準備
    private void space1(double x0, double y0, double x1, double y1) {
        X0 = x0; Y0 = y0;
        ratiox = WX / (x1 - x0);
        ratioy = WY / (y1 - y0);
    }
    
    // x-z ユーザー座標 (ワールド座標系) をウィンドウ座標 (デバイス座標系)
    private int wx(double x) {
      return (int)(ratiox * (x - X0));
    }
   
      //逆写像 
    private double wxinv(int ix) {
        return ix/ratiox + X0 ;
    }
    
    // x-z ユーザー座標 (ワールド座標系) をウィンドウ座標 (デバイス座標系)
    private int wy(double y) {
        return WY - (int)(ratioy * (y - Y0));
    }
   
      //逆写像 
    private double wyinv(int iy) {
        return (WY - iy)/ratioy + Y0 ;
    }
    
    // x-t,z-t 座標系の変換パラメーター
    private double ratioa,ratiob,A0,B0;
    private int WA = 350, WB =330;

    // x-t,z-t 描写範囲
    private double b_max =   1.5;
    private double b_min =  -0.5;

    private boolean IsIn2(double a, double b) {
        return (tt0 <= a && a <= tt && b_min <= b && b <= b_max);
    }

    // x-t,z-t 座標変換の準備
    private void space2(double b0,double b1) {
        A0 = tt0; B0 = b0;
        ratioa = WA / (tt - tt0);
        ratiob = WB / (b1 - b0);
    }
    
    // x-t,z-t ユーザー座標 (ワールド座標系) をウィンドウ座標 (デバイス座標系)
    private int wa(double a) {
      return (int)(ratioa * (a - A0));
    }
    
    // x-t,z-t ユーザー座標 (ワールド座標系) をウィンドウ座標 (デバイス座標系)
    private int wb(double b) {
        return WB - (int)(ratiob * (b - B0));
    }
    
    //数値の読み取り
    private void ReadFields(){
        q   = Double.valueOf(tf1.getText()).doubleValue();
        eps = Double.valueOf(tf2.getText()).doubleValue();
        xx0 = Double.valueOf(tf3.getText()).doubleValue();
        zz0 = Double.valueOf(tf4.getText()).doubleValue();
        cf  = Double.valueOf(tf5.getText()).doubleValue();
        tt  = Double.valueOf(tf6.getText()).doubleValue();
        x_min = Double.valueOf(tf8.getText()).doubleValue();
        x_max = Double.valueOf(tf9.getText()).doubleValue();
        y_min = Double.valueOf(tf10.getText()).doubleValue();
        y_max = Double.valueOf(tf11.getText()).doubleValue();
        x_margin = (x_max-x_min)/20;
        y_margin = (y_max-y_min)/20;
        maxn   = 1000;
        t = new double [maxn+1];
        x = new double [maxn+1];
        z = new double [maxn+1];
        x[0] = xx0;
        z[0] = zz0;
    }
  
    //関数の定義,ベクトル場の定義
    private double f(double x, double z)
    {
        return(x*(1-x)-cf*z*(x-q)/(q+x))/eps;
    }
    private double g(double x, double z)
    {

        return(x-z);
    }
 
   //ベクトルの設定
    
     //ベクトルの長さを均一にする
    double norm(double x, double y) {
        return Math.sqrt(x*x+y*y);
    }
    void arrow(double x0, double y0, double vx, double vy) {  /* ベクトルを引く
                                                                関数をつくる */
        double angle = 15 * Math.PI / 180.0;
        double cos = Math.cos(angle), sin = Math.sin(angle);
        double Norm = norm(vx, vy);
        double fx = 0.04 * (x_max - x_min), fy = 0.04 * (y_max - y_min);
        vx /= Norm; vy /= Norm;
        vx *= fx; vy *= fy;
        double xe = x0 + vx, ye = y0 + vy;
        vx *= - 0.2; vy *= -0.2;
        double vx1= vx*cos-vy*sin;     
        double vy1= vx*sin+vy*cos;
        double vx2= vx*cos+vy*sin;
        double vy2=-vx*sin+vy*cos;
        if (IsIn1(x0, y0) && IsIn1(xe, ye)) {
            line1(x0, y0, xe, ye);
            line1(xe, ye, xe+vx1, ye+vy1);
            line1(xe, ye, xe+vx2, ye+vy2);
        }
    }

   //線を引くメソッドを定義
    private void line1(double x0, double y0, double x1, double y1)
    {
        gra1.drawLine(wx(x0),wy(y0),wx(x1),wy(y1));
    }
    private void line2(double x0, double y0, double x1, double y1)
    {
        gra2.drawLine(wa(x0),wb(y0),wa(x1),wb(y1));
    }
     
    private double max(double a, double b) {
        if (a > b)
            return a;
        else
            return b;
    }
    private void report(String s) {
        String newStr = ta1.getText();
        if (newStr.length() < 1000)
            ta1.setText(newStr + s);
        else
            ta1.setText(s);
    }
    // RKF45 公式で t から t+h まで解を計算
    // 初期値 (x,z)
    // 結果 (next_x,next_z), 誤差見積もり *error
    private void rkf45_0(double t, double h, double x, double z) {
        double kx1,kx2,kx3,kx4,kx5,kx6,ky1,ky2,ky3,ky4,ky5,ky6;
        double X,Z,newx,newx2,newz,newz2;

        kx1=f(x,z);
        ky1=g(x,z);
        X = x+h*kx1/4.0; Z = z+h*ky1/4.0;
        kx2=f(X, Z);
        ky2=g(X, Z);
        X = x+h*(3.0*kx1+9.0*kx2)/32.0; Z = z+h*(3.0*ky1+9.0*ky2)/32.0;
        kx3=f(X,Z);
        ky3=g(X,Z);
        X = x+h*(1932.0*kx1-7200.0*kx2+7296.0*kx3)/2197.0;
        Z = z+h*(1932.0*ky1-7200.0*ky2+7296.0*ky3)/2197.0;
        kx4=f(X,Z);
        ky4=g(X,Z);
        X = x+h*(439.0/216.0*kx1-8.0*kx2+3680.0/513.0*kx3-845.0/4104.0*kx4);
        Z = z+h*(439.0/216.0*ky1-8.0*ky2+3680.0/513.0*ky3-845.0/4104.0*ky4);
        kx5=f(X,Z);
        ky5=g(X,Z);
        X = x+h*(-8.0/27.0*kx1+2.0*kx2-3544.0/2565.0*kx3+1859.0/4104*kx4-11.0/40.0*kx5);
        Z = z+h*(-8.0/27.0*ky1+2.0*ky2-3544.0/2565.0*ky3+1859.0/4104*ky4-11.0/40.0*ky5);
        kx6=f(X,Z);
        ky6=g(X,Z);
        newx=x+h*(16.0/135.0*kx1+6656.0/12825.0*kx3
                  +28561.0/56430.0*kx4-9.0/50.0*kx5+2.0/55.0*kx6);
        newz=z+h*(16.0/135.0*ky1+6656.0/12825.0*ky3
                  +28561.0/56430.0*ky4-9.0/50.0*ky5+2.0/55.0*ky6);
        newx2=x+h*(25.0/216.0*kx1+1408.0/2565.0*kx3
                   +2197.0/4104.0*kx4-1.0/5.0*kx5);
        newz2=z+h*(25.0/216.0*ky1+1408.0/2565.0*ky3
                   +2197.0/4104.0*ky4-1.0/5.0*ky5);
        next_x = newx;
        next_z = newz;
        error = norm(newx-newx2,newz-newz2);
    }
    // 時刻 t0 から tn まで解く
    // 初期値は (x0, z0)
    // 結果は xn, zn に入れて返す
    // 許容誤差限界を eps とする
    private void rkf45_1(double t0, double tn,
                         double x0, double z0,
                         double eps) {
        // 最大反復回数 maxiter
        int maxiter = 1000;
        // 計算機イプシロン程度の数
        double epsmac = 1e-15;
        int itend; // 最後のステップになるとき 1, そうでないとき 0
        double epsv = max(eps, epsmac); // まっとうな許容誤差限界
        double eb = Math.abs(epsv / (tn - t0)); // 単位時間あたりの許容誤差限界
        Format f_f = new Format("%6.3f");
        Format f_e = new Format("%6.3e");
        Format f_g = new Format("%6.3g");
        double t, h;

        //
        t = t0;
        h = tn - t0;
        itend = 0;

        ta1.setText(""+f_f.form(t0)+"から"+f_f.form(tn)+"まで解く\n"
                   +"(x,z)=("+f_f.form(x0)+","+f_g.form(z0)+")\n");
        for (int iter=1; iter<=maxiter; iter++) {
            rkf45_0(t, h, x0, z0);
            xn = next_x; zn = next_z;
            if (verbose)
                report("rkf45_0()から戻る\n iter="+iter
                       +",(x,z)=("+f_g.form(xn)+","+f_g.form(zn)+"),\n"
                       +"error="+f_e.form(error)+"\n");
            if (error <= Math.abs(eb * h)) { // 計算がちゃんとできた
                if (itend == 1) { // 最後の計算だったら
                    report(" ちゃんと最後まで行った。\n");
                    return;
                }
                itend = 0;
                t += h;
                h *= 0.9 * Math.sqrt(Math.sqrt(Math.abs(eb * h) / error));
                if (verbose)
                    report("t="+f_g.form(t)+"\n"+" 新しいh="+f_e.form(h)+"\n");
                if (Math.abs(h) > Math.abs(tn-t)) {
                    // もう残り時間は短くて、|h| 未満
                    h = tn - t;
                    itend = 1;
                }
                x0 = xn;
                z0 = zn;
            }
            else {
                // 推定誤差時間刻み幅 h あたりの許容誤差よりも大きい場合
                if (error > epsv) {
                    h *= 0.1;
                    report("h が大きいので小さくした。"+"\n"
                           +"新 h="+f_e.form(h)+"\n");
                }
                else {
                    h *= 0.9 * Math.sqrt(Math.sqrt(Math.abs(eb * h) / error));
                    if (verbose)
                        report("新しいh="+f_e.form(h)+"\n");
                    if (Math.abs(h) < epsmac) {
                        // 非常に小さい時間刻み幅が出てきた。もう無理だ！
                        // 警告を発すべきである mk
                        report(" h="+f_e.form(h)+" は小さくなりすぎ！\n");
                        return;
                    }
                }
            }
        } // for ループの終わり
        report("rkf45_1(): 最後まで計算したけれど...\n");
    } // rkf45_1() の終わり

    private void rkf45(double ts, double te,
                       double []t, double []x, double []z,
                       double eps) {
        int i;
        n = 1000;
        if (n > maxn)
            n = maxn;
        double dt = (te - ts) / n;
        for (i = 0; i <= n; i++)
            t[i] = ts + i * dt;
        for (i = 0; i < n; i++) {
            rkf45_1(t[i], t[i+1], x[i], z[i], eps / n);
            x[i+1] = xn;
            z[i+1] = zn;
        }
    }
   
    //パラメーター表示   
    public void init() {
        setCursor(new Cursor(Cursor.HAND_CURSOR));

        //tf la bt ta のレイアウト
        setLayout(null);
        add(la1 = new Label("ｑの値:"));
        la1.setBounds (650,10,90,20);
        add(la2 = new Label("εの値:"));
        la2.setBounds (835,10,90,20);
        add(la3 = new Label("初期値のｘ座標:"));
        la3.setBounds (650,30,90,20);
        add(la4 = new Label("初期値のｚ座標:"));
        la4.setBounds (835,30,90,20);
        add(la5 = new Label("ｆの値:"));
        la5.setBounds (650,50,90,20);
        add(la6 = new Label("ｔ座標の上限:"));
        la6.setBounds (835,50,90,20);
        //add(la7 = new Label("ｆの値:"));
        //la7.setBounds (650,70,90,20);
//      add(la8 = new Label("区間の分割数:"));
        //      la8.setBounds (835,70,90,20);
        add(la7 = new Label("平衡点の座標:"));
        la7.setBounds (650,130,80,20);
        add(la8 = new Label(" x の範囲"));
        la8.setBounds (650,155,50,20);
        add(la9 = new Label("〜"));
        la9.setBounds (750,155,20,20);
        add(la10 = new Label(" z の範囲"));
        la10.setBounds (650,175,50,20);
        add(la11 = new Label("〜"));
        la11.setBounds (750,175,20,20);
        add(ta1 =new TextArea(5,7));        //TextArea
        ta1.setBounds (650,200,175,100);
            
        //初期値などの代入
        add(tf1 = new TextField("" +8E-4)); 
        tf1.setBounds (745,10,80,20);
        add(tf2 = new TextField("" +1E-2)); 
        tf2.setBounds (930,10,80,20);
        add(tf3 = new TextField("" +0.01)); 
        tf3.setBounds (745,30,80,20);
        add(tf4 = new TextField("" +0.01)); 
        tf4.setBounds (930,30,80,20);
        add(tf5 = new TextField("" +0.6)); 
        tf5.setBounds (745,50,80,20);
        add(tf6 = new TextField("" +init_tt)); // tt
        tf6.setBounds (930,50,80,20);
        //add(tf7 = new TextField("" +0.6));
        //tf7.setBounds (745,70,80,20);
        //add(tf8 = new TextField("" +100000));
        //tf8.setBounds (930,70,80,20);
        add(tf7 = new TextField(""));
        tf7.setBounds (730,130,350,20);
        add(tf8 = new TextField("" +0.0));
        tf8.setBounds (700,155,50,20);
        add(tf9 = new TextField("" +1.0));
        tf9.setBounds (770,155,50,20);
        add(tf10 = new TextField("" +0.0));
        tf10.setBounds (700,175,50,20);
        add(tf11 = new TextField("" +1.0));
        tf11.setBounds (770,175,50,20);

        //ボタンの設定

          //スクリーン 1 Startボタン
        add(bt1  = new Button("Start 1"));
        bt1.addActionListener(this);
        bt1.setBounds (650,70,175,20);

          //ベクトル場 start ボタン
        add(bt2  = new Button("ベクトル場表示"));
        bt2.addActionListener(this);
        bt2.setBounds (650,90,175,20);
        
          //スクリーン 1 Clear ボタン
        add(bt3  = new Button("Clear 1"));
        bt3.addActionListener(this);
        bt3.setBounds (650,110,175,20);

          //スクリーン 2 Startボタン
        add(bt4  = new Button("Start 2"));
        bt4.addActionListener(this);
        bt4.setBounds (825,70,175,20);
        
          //スクリーン 2 Clearボタン
        add(bt5  = new Button("Clear 2"));
        bt5.addActionListener(this);
        bt5.setBounds (825,90,175,20);
    
        //スクリーンの設定

           //スクリーン 1 の確保
        img1=createImage(WX,WY);
        gra1=img1.getGraphics();
        gra1.setColor(new Color(230,255,230));
        gra1.fillRect(0,0,WX,WY);
        
           // 座標変換
        ReadFields();
        space1(x_min-x_margin, y_min-y_margin,
               x_max+x_margin, y_max+y_margin);
        
           // ラベル
        gra1.setColor(Color.black);
        gra1.drawString("オレゴネータ（２変数版 ）の解曲線",200,30);
        
           //スクリーン 2 の確保
        img2=createImage(WA,WB);
        gra2=img2.getGraphics();
        gra2.setColor(new Color(230,255,230));
        gra2.fillRect(0,0,WA,WB);
   
           //ラベル
        gra2.setColor(Color.red);
        gra2.drawString("ｔ−ｘ",100,10);
        gra2.setColor(Color.orange);
        gra2.drawString("ｔ−ｚ",150,10);
              
        gra2.setColor(Color.black);
        gra2.drawString("のグラフ",200,10);

        
        addMouseListener(this);
    
    }
        
    public void paint(Graphics g){
        update(g);
    }
    public void update(Graphics g){
        g.drawImage(img1,0,0,this);
        g.drawImage(img2,655,320,this);
    }

    //ボタン操作 
    public void actionPerformed(ActionEvent e) {

        //ボタン 1
        if(e.getSource() == bt1) {          
            // 初期値などを取得
            ReadFields();

            // RKF 45 で解を計算
            rkf45(0.0, tt, t, x, z, eps_tol);

            // 解曲線を描く
            space1(x_min-x_margin, y_min-y_margin,
                   x_max+x_margin, y_max+y_margin);
            gra1.setColor(Color.blue);
            for(int j=1; j<=n; j++) {
                if (IsIn1(x[j-1], z[j-1]) && IsIn1(x[j], z[j])) {
                    line1( x[j-1], z[j-1], x[j], z[j]);
                }
            }
         
            gra1.setColor(Color.lightGray);
            for (double i=x_min; i<=x_max ; i=i+0.1){
                line1( i, y_min, i, y_max);
            }
            for (double i=y_min; i<=y_max; i=i+0.1){
                line1( x_min, i, x_max, i);
            }

            // ラベルを描く
            gra1.setColor(Color.black);
            line1( x_min, 0.0, x_max, 0.0);
            line1( 0.0, y_min, 0.0, y_max);
            gra1.drawString("0",wx(0.0),wy(0.0));
            gra1.drawString("オレゴネータ（２変数版 ）の解曲線",200,30);
            for (double i=0.5; i<=x_max; i=i+0.5){
                gra1.drawString(""+i,wx(i),wy(0.0));
            }
            for(double i=0.5; i<=y_max; i=i+0.5){
                gra1.drawString(""+i,wx(0.0),wy(i));
            }
            repaint();

            //平衡点を計算する
            ReadFields();
            double x_h =(1-cf-q+Math.sqrt(Math.pow(cf+q-1,2)+4*q*(1+cf)))/2;

            //平衡点をプロット
            gra1.setColor(Color.red);
            double h1= x_h-(x_max-x_min)/100;
            double h2= x_h+(x_max-x_min)/100;
            double h3= x_h-(x_max-x_min)/100;
            double h4= x_h+(x_max-x_min)/100;
            line1( h1, x_h, h2, x_h);
            line1( x_h, h3, x_h, h4);

            //平衡点を表示させる
            String he=Double.toString(x_h);
            tf7.setText("("+ he + "," + he +")");
          
        } // end if

        //ボタン 2
          //ベクトル場表示
        else if(e.getSource() == bt2){
          ReadFields();
          gra1.setColor(Color.green);
          double dx = (x_max - x_min) / 20, dy = (y_max - y_min) / 20;
          for(double x=x_min; x<=x_max; x += dx){   /* 矢印を書く */
              for(double y=y_min; y<=y_max; y += dy){  
                  arrow(x, y, f(x,y), g(x,y));
              }
          }
          repaint();
        }
        
        //ボタン 3  
          // 相曲線とベクトル場をクリア
          // ...相図をクリアしてグリッドなど描き直し
        else if(e.getSource() == bt3) {
            gra1.setColor(new Color(230,255,230));
            gra1.fillRect(0,0,WX,WY); 
        
            // ラベルを描く
            gra1.setColor(Color.black);
            gra1.drawString("オレゴネータ（２変数版 ）の解曲線",200,30);

            repaint();

            //初期値をはじめの値にする
            tf3.setText("" +0.01);
            tf4.setText("" +0.01);
        
            //平衡点の値を消す
            tf7.setText("" );

            ta1.setText("");
        }
        
        //ボタン 4
           // 解曲線 (t,x(t)), (t,z(t)) を描く
        else if(e.getSource() == bt4) { 
            double TT,TT0;
            // 初期値などを取得
            ReadFields();
            space2(b_min, b_max);
                      
            //座標軸
            gra2.setColor(Color.black);
            line2( tt0, 0.0, tt, 0.0);
            line2( tt0, b_min, tt0, b_max);
            gra2.drawString( ""+tt0 ,wa(tt0),wb(0.0));
            gra2.drawString( "1.0",wa(tt0),wb(1.0));
            gra2.drawString(""+tt,wa(0.9*tt),wb(0.0));
            gra2.drawString(" x,z",wa(tt0),wb(1.4));

            // RKF45 で解を計算
            rkf45(0.0, tt, t, x, z, eps_tol);

            for (int j=1; j<=n; j++) { 
                  if (IsIn2(t[j-1], x[j-1]) && IsIn2(t[j], x[j])) {
                      gra2.setColor(Color.red);
                      line2(t[j-1],x[j-1],t[j], x[j]);
                  }
                  if (IsIn2(t[j-1], z[j-1]) && IsIn2(t[j], z[j])) {
                      gra2.setColor(Color.orange);
                      line2(t[j-1], z[j-1], t[j], z[j]);
                  }
                  if (Math.abs(t[j]) >= Math.abs(tt))
                      break;
                  
            }
            repaint();
        }  // end if

        //ボタン 5      
           // 解曲線 (t,x(t)), (t,z(t)) 消去
        else if(e.getSource() == bt5) {
            gra2.setColor(new Color(230,255,230));
            gra2.fillRect(0,0,WA,WB);
           
         gra2.setColor(Color.red);
              gra2.drawString("ｔ−ｘ",100,10);
              gra2.setColor(Color.orange);
              gra2.drawString("ｔ−ｚ",150,10);

              gra2.setColor(Color.black);
              gra2.drawString("のグラフ",200,10);
              
              repaint();
        }
        
    }
    
 
    public void mousePressed(MouseEvent e){}
    public void mouseReleased(MouseEvent e){} 
    public void mouseClicked(MouseEvent e){

        //マウスによる初期値入力
        // if(e.getModifiers()==){
            space1(x_min-x_margin,y_min-y_margin,x_max+x_margin,y_max+y_margin);
            int xx=e.getX();
            int zz=e.getY();
            double x_m=wxinv(xx);
            double z_m=wyinv(zz);
            String xx0=Double.toString(x_m);
            String zz0=Double.toString(z_m);
            tf3.setText(xx0);
            tf4.setText(zz0);

              // 初期値などを取得
            ReadFields();
            // RKF45 で解を計算
            rkf45(0.0, tt, t, x, z, eps_tol);
            // 解曲線を描く
            gra1.setColor(new Color(0,0,200));
            for (int j=1; j<=n; j++) {
                if (IsIn1(x[j-1], z[j-1]) && IsIn1(x[j], z[j])) {
                    line1( x[j-1], z[j-1], x[j], z[j]);
                }
            }
           
            repaint();
    }  // end MouseClick
    public void mouseEntered(MouseEvent e){}
    public void mouseExited(MouseEvent e){}

}