/* metacircle morphing by using Bezier Clipping by T.N '97 11/17
  modified on 97' 11/15, morph(11/20), mesh(11/23)
                            */
import java.applet.Applet;
import java.awt.*;
import java.io.*;
import java.net.URL;
import java.util.StringTokenizer;

/* *********************************************************: */
public class circmor extends Applet {
	Point[] p = new Point[10];
	Point[] p1 = new Point[10];
	Point PntC;
	int n = 0;
	int Nc=4;
	int ndiv = 100;
	int disp_con = 1;
	int disp_mod = 0;
	int curve1=0;
	double Rmax=300.;
	double TOLEL= 0.005;
	int multi = 0;
	int Ncurve=-1;
    double tintt;
    int nsplit=0,n_clip=0,iproces;
	Button b1,b2,b3,b4,b5;
	Choice cmod = new Choice();
	Choice cn1 = new Choice(); //cn1に項目の追加
	Choice cn2 = new Choice();
	Panel p3 = new Panel();
	Label l1,l6;
	TextField t1;
	Point prev = new Point(0, 0);
	int idrag=0;
	double select_u=-1.;

	Curve3[] Bez=new Curve3[10];
	Curve3[] Bez_back=new Curve3[10];
	Curve3[] Bez1=new Curve3[10];
	Curve3[] Bez_s=new Curve3[10];
	Curve3[] Bez_d=new Curve3[10];

    int closestCurv=0;
	double DISmin=100000.;
	Point closestPnt = new Point(0,0);

 	int x_origin,y_origin;
	double	scale=1.;
    int  view_chang=1;
     int findp=0;

	 int itest =1;
	 int MAXdeg = 20;
	 double[][] binom = new double[MAXdeg][MAXdeg];
	 int idebug=0;int nbar=0; 

	 int addcirc=0;
	 int ncirc=0,ncirc_s=0,ncirc_d=0;
	 Metaball metaball[];
	 Metaball metaball_s[];
	 Metaball metaball_d[];

	 Vec2 circp[][];
	 Vec2 circp_s[][];
	 Vec2 circp_d[][];

	 int inc_circ=0;
	 String datafile=" ";
	 double xcenter= 0., ycenter = 0.;
	 double half1 = 0.7;
	 int fill_crc=0;
     int disp_anime=0;
	 int morph=0;
	 int cores[];
	 int Nmesh=0,Mmesh=0,mesh_move=0;

/* **************************************************************: */
	public void init() {
	String str;
	Dimension appsize;

		setBackground(new Color( 140, 140,150));

		setForeground(Color.red);
		appsize = size();

		setLayout(new BorderLayout(0,0));
		CheckboxGroup gc = new CheckboxGroup();
        Checkbox c2= new Checkbox("add",gc, true);
        Checkbox c1= new Checkbox("deform",gc, false);

		Panel p = new Panel();
        p.setBackground(new Color(40,40,125));

		p.add(b1 = new Button("New Curve"));
		p.add(b2 = new Button("Add Curve"));
		p.add(new Label("Mode;"));
		p.add(c1);p.add(c2);
		p.add(b3 = new Button("dbg"));
        add("South",p);

		//cmodに項目の追加
		cmod.addItem("Input Mode");
		cmod.addItem("HourceData");
		cmod.addItem("TigerData");
		cmod.addItem("CircleFill");
		cmod.addItem("Clear");
		cmod.addItem("Animation");
		cmod.addItem("Data Print");
		cmod.addItem("Mesh");

//		Choice cn1 = new Choice(); //cn1に項目の追加
		cn1.addItem("degree 3 Bezier Curve");
		cn1.addItem("degree 4 Bezier Curve");
		cn1.addItem("degree 5 Bezier Curve");
		cn1.addItem("degree 6 Bezier Curve");
		cn1.addItem("degree 7 Bezier Curve");
		cn1.addItem("degree 8 Bezier Curve");

		//p3に項目の追加
        p3.setBackground(new Color(40,40,125));
		p3.add(cmod);
		p3.add(cn1);

	//	p3.add(l6 = new Label("R="));
	//	p3.add(t1 = new TextField("100",3));
		p3.add(b4 = new Button("back"));
		p3.add(b5 = new Button("morph"));

        Checkbox d2= new Checkbox("test");
		p3.add(d2);

        add("North",p3);
		metaball = new Metaball[500];
		circp = new Vec2[500][5];

	if ((str = getParameter("disp_mod")) != null) {
	    disp_mod= Integer.parseInt(str) ; }
 		resize(350,340);
		showStatus("Click Draw bottun!");
	    x_origin = appsize.width/2; y_origin = appsize.height/2;
		binit(18);
    }
/* *************** 描画関数 *****************************: */
    public void paint(Graphics g) {
	 int i,ix; //,nint=0
 //    int ndeg1=4;

	 if(idrag==0) {
		g.setColor(Color.yellow);
		 for(i = 0; i < n; i++) g.fillOval(p[i].x-3, p[i].y-3, 7, 7);
 		 if(n>1) {g.setColor(Color.black);
		 for(i = 1; i < n; i++) g.drawLine(p[i-1].x,p[i-1].y,p[i].x,p[i].y); }
	 }
	 if(Nmesh>0) meshDisp(g);
	 if(ncirc>0) { circDisp(g); } //円の描画		 
	 if(morph >= 2) {preMorp(g); return;} // モーフイング

	 if(Ncurve>0 ) {
		if(idrag==0 || findp==0) {
		   for(int k=0;k<Ncurve;k++) {
			 for(i = 0; i < Nc; i++) Bez[k].p3[i] = tranfrm(Bez[k].P3[i]);
			 Bez[k].Polygon3(Nc,g); 
			 g.setColor(Color.red); Bez[k].draw3(Nc,g);
//		    System.err.println(" k< Nc="+Nc+" k="+k+" Nc="+Ncurve);
		   }
	   }
	  }

	   if(curve1==0) {
	     if(n==Nc) {	//g.setColor(Color.red);
			 if( Ncurve<0)  Ncurve=0;
		    for(i=0;i<Nc;i++) p1[i]=p[i];
			Bez[Ncurve] = new Curve3();
		    Bez[Ncurve].SetOrg(x_origin,y_origin);
	        for(i=0;i<Nc;i++) Bez[Ncurve].p[i]= new Point(p[i].x,p[i].y);
			Bez[Ncurve].copy(Nc); 
		    for(i=0; i<Nc; i++) Bez[Ncurve].p3[i] = tranfrm(Bez[Ncurve].P3[i]);
		//	if(idrag == 1) Bez[Ncurve].Polygon3(Nc,g);
			 g.setColor(Color.red); Bez[Ncurve].draw3(Nc,g);
//			System.err.println(" Nc="+Nc+" p1="+p1[0].x+" Nc="+Ncurve);

			Bez_back[Ncurve] = new Curve3();
		    Bez_back[Ncurve].SetOrg(x_origin,y_origin);
	        for(i=0;i<Nc;i++) Bez_back[Ncurve].p[i]= new Point(p[i].x,p[i].y);
			Bez_back[Ncurve].copy(Nc); 

			n=0; curve1=1;
          }
       }         // end :if(n=NC) 
       else {   // for curve=1,2
		 // plot point
		 if(n==1) { g.setColor(Color.yellow);
			if(idrag == 1 && disp_anime==0) {
             for(int k=0;k<=Ncurve; k++) {
			   for(i=0; i<Nc; i++) Bez[k].p3[i] = tranfrm(Bez[k].P3[i]);
//             曲線表示
               Bez[k].Polygon3(Nc,g); g.setColor(Color.red); Bez[k].draw3(Nc,g);
			 }
   //         System.err.println("drag Nc="+Nc+"Ncuve="+Ncurve);
			}
			else if(Ncurve>0) {
              for(int k=0;k<Ncurve; k++) {
			   for(i=0; i<Nc; i++) Bez[k].p3[i] = tranfrm(Bez[k].P3[i]);
                Bez[k].Polygon3(Nc,g);
				g.setColor(Color.red); Bez[k].draw3(Nc,g);
			   }
			  }
              if(findp==1) curve1 = 2;	
		 }
		 // plot degree Nc curve 
		       for (int k=0;k<=Ncurve;k++){
			     for(i=0; i<Nc; i++) Bez[k].p3[i] = tranfrm(Bez[k].P3[i]);
				 if(disp_anime==0) {
			     g.setColor(Color.gray);Bez[k].Polygon3(Nc,g);
				 g.setColor(Color.green); Bez[k].draw3(Nc,g);}
			   }
			   if(idrag==1) { // ドラッグ中の円描画
				  g.setColor(Color.black);  defrmcirc(g); }
			   n=0;
    //        System.err.println("Drag Nc="+Nc+"Ncuve="+Ncurve+"drag="+idrag);
	   }
	   if(disp_mod == 2) {  // データ表示
			 for(int k=0;k<=Ncurve;k++) DispData(k,g); return; }
	   if(findp==1 && curve1==2) { CPfind(g); } // 最近点の抽出
   }
 /* ****************** Display Circles ***********************: */
	public void circDisp(Graphics g) {
		int i;
  
	  g.setColor(Color.red);if(idrag==1) g.setColor(Color.gray);
	  if(idrag==0) { if(fill_crc==1) g.setColor(Color.black);
		for (i=1; i<=ncirc; i++) {
				int ir=(int)(metaball[i].r*half1);

				if(fill_crc==1)	g.fillOval((int)(metaball[i].cen.x)-ir,
					(int)(metaball[i].cen.y)-ir,ir+ir,ir+ir);
				else
				g.drawOval((int)(metaball[i].cen.x)-ir,
					(int)(metaball[i].cen.y)-ir,ir+ir,ir+ir);
		}
	}
	if(inc_circ==0) {
		    g.setColor(Color.white); g.drawString("Waite a moment:",20,65);
			for (i=1; i<=ncirc; i++) {
                Point pc = new Point((int)(metaball[i].cen.x),(int)(metaball[i].cen.y));
				circpos(i,pc,g);
			}
         inc_circ=1;
         g.drawString("Please deform",20,75);
	}
}
 /* ****************** Mesh ***********************: */
	public void meshDisp(Graphics g) {
	int i,j;
	double x1,x2,y1,y2;
	double Width =270.,Hight=180.;
	int Yoff=70,Xoff=50;
	
	if(findp==0) { g.setColor(Color.blue);
	  for(i=0;i<=Nmesh;i++){
		for(j=0;j<Mmesh; j++) {
			x1 = Width*(double)i/(double)Nmesh+Xoff;
			y1 = Hight*(double)j/(double)Mmesh+Yoff;
			y2 = Hight*(double)(j+1)/(double)Mmesh+Yoff;
			g.drawLine((int)x1, (int)y1, (int)x1,  (int)y2);
		}
	  }
	  for(j=0;j<= Mmesh;j++){
		for(i=0;i<Nmesh; i++) {
			x1 = Width*(double)i/(double)Nmesh+Xoff;
			x2 = Width*(double)(i+1)/(double)Nmesh+Xoff;
			y1 = Hight*(double)j/(double)Mmesh+Yoff;
			g.drawLine((int)x1, (int)y1, (int)x2, (int)y1);
		}
	  }
	}
	else {
		if(mesh_move == 0) {
			int nsize= (Nmesh+1)*(Mmesh+1)+1;
		circp_s = new Vec2[nsize][1];
		int ii =0;
	  for(j=0;j<= Mmesh;j++){
		for(i=0;i<=Nmesh; i++) {			
			x1 = Width*(double)i/(double)Nmesh+Xoff;
			y1 = Hight*(double)j/(double)Mmesh+Yoff;
           Point pc = new Point((int)(x1),(int)(y1));
		   meshPos(ii,pc,g); ii++;

		}
	  }
		        mesh_move=1; 
		}
	 else {
		 g.setColor(Color.red);
		int ii =0;
	   for(j=0;j<= Mmesh;j++){
		 for(i=0;i<=Nmesh; i++) {
 		   Point2 ps = new_point(ii,circp);
		   circp_s[ii][0] = new Vec2(ps.x, ps.y);		 
		   g.fillOval((int)ps.x+x_origin-2, (int)ps.y+y_origin-2, 5, 5);
			ii++;
		 }
	   }
	   g.drawString("mesh points:"+ii,15,45);
	 }
	}
}
	/* ***********************************************************: */
	public void meshPos(int i_mesh, Point p,Graphics g){
	Point3[] P3 = new Point3[10];
	 Point[] p2 = new Point[10];
	 Point p0;
	 Vec2 u_dist= new Vec2(0.,0.);
	 Point3 P0 =new Point3(0.,0.,0.);
	 Point3 P1 =new Point3(0.,0.,0.);
	  int i;
	  double a,b,c,d,f;

       int ndeg1 = Nc; 
	  for(int k=0;k<=Ncurve;k++){ 
				  DISmin = 100000.;
  		        for(i=0;i<Nc;i++) p2[i]=Bez[k].p[i];
			    for(i=0;i<Nc;i++) P3[i] = new Point3((int)(Bez[k].p3[i].x)+x_origin,
				 (int)(Bez[k].p3[i].y)+y_origin,Bez[k].p3[i].z);
				u_dist = new Vec2(0.,0.);
                p0=clospnt3(ndeg1-1,P3,p2, p, k, DISmin,u_dist,g); //B.C.法利用

				  circp[i_mesh][k] = new Vec2(u_dist.u, u_dist.v);
				  double u = circp[i_mesh][k].u;
				  P0 = (Point3)Bez[k].eval(Nc,u);  
				  P1 = (Point3)Bez[k].eval(Nc,u+0.01);
				  a = -(P1.z-P0.z); b =(P1.y-P0.y);
				  d = Math.sqrt(a*a+b*b);
				  a = a/d; b=b/d; c = -(a*(double)p0.x +b*(double)p0.y);
                  f = a*(double)p.x +b*(double)p.y+c;
				  if(f<0.) circp[i_mesh][k].v = -u_dist.v;
/*				  double x = P0.y+a* circp[i_mesh][k].v;
				  double y = P0.z+b* circp[i_mesh][k].v;
*/
//				g.setColor(Color.pink); g.fillOval((int)P0.y+x_origin-2, (int)P0.z+y_origin-2, 3, 3);
	 }
   }

/* ****************** 最近点の抽出 ***********************: */
	public void CPfind(Graphics g) {
	int i;
    int ndeg1=4;
	 Point[] p2 = new Point[10];
	 Point3[] P3 = new Point3[10];
	 Vec2 u_dist = new Vec2(0.,0.);

			ndeg1 = Nc; DISmin = 100000.; select_u=-1.;
	          		 // 最後の曲線をまず処理
		   for(i=0;i<Nc;i++) P3[i] = new Point3((int)(Bez[Ncurve].p3[i].x)+x_origin,
				      (int)(Bez[Ncurve].p3[i].y)+y_origin,Bez[Ncurve].p3[i].z);
              p[2]=clospnt3(ndeg1-1,P3,p1, p[0], Ncurve, Rmax, u_dist,g); //B.C.法利用
		   g.setColor(Color.white); g.fillOval(p[2].x-3, p[2].y-3, 7, 7);
		   if(Ncurve>0) {
			  for(int k=0;k<Ncurve;k++){ // 他の曲線
  		        for(i=0;i<Nc;i++) p2[i]=Bez[k].p[i];
			    for(i=0;i<Nc;i++) P3[i] = new Point3((int)(Bez[k].p3[i].x)+x_origin,
				 (int)(Bez[k].p3[i].y)+y_origin,Bez[k].p3[i].z);
                p[2]=clospnt3(ndeg1-1,P3,p2, p[0], k, DISmin,u_dist,g); //B.C.法利用
			    g.setColor(Color.white); g.fillOval(p[2].x-3, p[2].y-3, 7, 7);
		        }
		   }
			  if(disp_anime==0) {g.setColor(Color.white); 
			  Bez[closestCurv].draw3(Nc,g); } // 選択された曲線再表示
              curve1=1;  n=0;
}
/* ******************デフォーメーション *********************: */
	public void defrmcirc(Graphics g) {
	int i;

	 g.setColor(Color.red);
	 if(fill_crc == 1) g.setColor(Color.black);
	 for (i=1; i<=ncirc; i++) {				
				Point2 ps = new_point(i,circp);
				int ir=(int)(metaball[i].r*half1);
				if(fill_crc == 1) 
					g.fillOval((int)ps.x+x_origin-ir,(int)ps.y+y_origin-ir,ir*2,ir*2);
				else g.drawOval((int)ps.x+x_origin-ir,(int)ps.y+y_origin-ir,ir+ir,ir+ir);
	}
}
/* ********* データロード FileInp **********************: */
	public void FileInp() {
		int i;
		try {	
			InputStream is = new URL(getDocumentBase(), datafile).openStream();
			DataInputStream dis = new DataInputStream(is);
//			BufferedReader dis = new BufferedReader (new InputStreamReader(is));
			try {
				StringTokenizer st;
				double x, y, z,r,v;
				String comment;

				comment = dis.readLine();
				st = new StringTokenizer(dis.readLine());
//				n = (int)Double.valueOf(st.nextToken()).doubleValue();
				ncirc = Integer.valueOf(st.nextToken()).intValue();
				System.err.println("n=" + ncirc); 
	//			metaball = new Metaball[500];
	//			metaball = new Metaball[ncirc+1];

				for (i=1 ; i<=ncirc ;i++) {
					st = new StringTokenizer(dis.readLine());
					x = scale*Double.valueOf(st.nextToken()).doubleValue();
					y = scale*Double.valueOf(st.nextToken()).doubleValue();
					z = scale*Double.valueOf(st.nextToken()).doubleValue();
					r = scale*Double.valueOf(st.nextToken()).doubleValue();
					v = Double.valueOf(st.nextToken()).doubleValue();
					metaball[i] = new Metaball(x+xcenter,y+ycenter,z,r,v);
				}
			} 
			catch(IOException e) {
				//gr.drawString("Data read err:=",20,110);
				System.err.println("Can't access datafile " + datafile);
			}
			finally {
				is.close();
			}
		}
		catch(IOException e) {
			System.err.println("Closing error");
			System.exit(1);
		}
	}
/* *********************************************************: */
		public Point3 tranfrm(Point3 point) {
		return new Point3(point.y,point.z,1.);
	}
/* ***********************************************************: */
	public void circpos(int i_ball, Point p,Graphics g){
	Point3[] P3 = new Point3[10];
	 Point[] p2 = new Point[10];
	 Point p0;
	 Vec2 u_dist= new Vec2(0.,0.);
	 Point3 P0 =new Point3(0.,0.,0.);
	 Point3 P1 =new Point3(0.,0.,0.);
	  int i;
	  double a,b,c,d,f;

       int ndeg1 = Nc; 
	  for(int k=0;k<=Ncurve;k++){ 
				  DISmin = 100000.;
  		        for(i=0;i<Nc;i++) p2[i]=Bez[k].p[i];
			    for(i=0;i<Nc;i++) P3[i] = new Point3((int)(Bez[k].p3[i].x)+x_origin,
				 (int)(Bez[k].p3[i].y)+y_origin,Bez[k].p3[i].z);
				u_dist = new Vec2(0.,0.);
                p0=clospnt3(ndeg1-1,P3,p2, p, k, DISmin,u_dist,g); //B.C.法利用

				  circp[i_ball][k] = new Vec2(u_dist.u, u_dist.v);
				  double u = circp[i_ball][k].u;
				  P0 = (Point3)Bez[k].eval(Nc,u);  
				  P1 = (Point3)Bez[k].eval(Nc,u+0.01);
				  a = -(P1.z-P0.z); b =(P1.y-P0.y);
				  d = Math.sqrt(a*a+b*b);
				  a = a/d; b=b/d; c = -(a*(double)p0.x +b*(double)p0.y);
                  f = a*(double)p.x +b*(double)p.y+c;
				  if(f<0.) circp[i_ball][k].v = -u_dist.v;
				  double x = P0.y+a* circp[i_ball][k].v;
				  double y = P0.z+b* circp[i_ball][k].v;

//				g.setColor(Color.pink); g.fillOval((int)P0.y+x_origin-2, (int)P0.z+y_origin-2, 3, 3);

				  if(i_ball%10==0 || u_dist.v>120.) {
//				  System.err.println("k="+k+" ball="+i_ball+" u="+circp[i_ball][k].u+" r="+circp[i_ball][k].v+" f="+f);
//          System.err.println("ball="+i_ball+" u="+circp[i_ball][k].u+" x="+P0.y+" y="+P0.z);
//				  System.err.println("x="+x+" y="+y);
//          g.drawLine((int)x+x_origin, (int)y+y_origin, (int)P0.y+x_origin,  (int)P0.z+y_origin);
				  }
//				  if(k==1) {g.setColor(Color.gray);
//				g.fillOval((int)x+x_origin-1, (int)y+y_origin-1, 2, 2);}
//				else {g.setColor(Color.green); 
//				g.fillOval((int)x+x_origin-1, (int)y+y_origin-1, 4, 4);}
/*

System.err.println("p1x="+p2[0].x+"p1y="+p2[0].y);
System.err.println("P1x="+P3[0].x+"P1y="+P3[0].y);
*/
  //     System.err.println("xorg="+x_origin+" yorg="+y_origin);
    //     System.err.println("a="+a+" b="+b+" c="+c+" d="+d+" f="+f);
	 }
   }
/* **************************************************************: */
	public void DispData(int k,Graphics gr) {
		gr.setColor(Color.black);
     for(int i=0;i<Nc;i++) {
		 int kk = k*80+80;
		gr.drawString("i = " + i, 19, 80+i*14); 
		gr.drawString("x = " + Bez[k].P3[i].z, 50, kk+i*14); 
		gr.drawString("y = " + Bez[k].P3[i].y, 120, kk+i*14); 
		gr.drawString("z = " + Bez[k].P3[i].x, 190, kk+i*14); }
    }
/* ************ Ｂezier clipping による最近点検出 ********: */
	public Point clospnt3(int ndeg,Point3 Pnt[],Point pnt[], Point C, 
		int Cnum, double dRmax, Vec2 u_dist, Graphics g)   {
		Point  PntQ= new Point(0,0);
		Point[] pnt1 = new Point[11];	Point[] pnt2 = new Point[10];
		Point3[] Pnt1 = new Point3[11];  Point3[] Pnt2 = new Point3[10];

		double r2,u_min=0.;  int i,nint;
        double[] t1= new double[2];
        double[] x1= new double[11]; double[] y1= new double[11];
        double[] x2= new double[10]; double[] y2= new double[10];
        double[] w1= new double[11]; double[] w2= new double[10]; 

        double[] xint= new double[4];
        double[] xint1= new double[4];

		for( i=0;  i <=ndeg; i++) { // 2次元 平行移動
		pnt1[i] =  new Point(pnt[i].x-C.x,pnt[i].y-C.y);}
		for( i=0;  i <=ndeg; i++) { // ３次元 平行移動
		Pnt1[i] =  new Point3(Pnt[i].x-(double)C.x,Pnt[i].y-(double)C.y,Pnt[i].z);}

		PntC = new Point(C.x,C.y);
        if(disp_con==0){ // debug時のみ
	 	 for( i=0;  i <ndeg; i++) { //ホドグラフ計算
		 pnt2[i] =  new Point(pnt[i+1].x-pnt[i].x,pnt[i+1].y-pnt[i].y); }
		 for(i=1;  i <ndeg; i++) {
          g.drawLine(pnt2[i-1].y+C.x,-pnt2[i-1].x+C.y,pnt2[i].y+C.x,-pnt2[i].x+C.y);}
          g.setColor(Color.darkGray);
          g.drawLine(C.x, C.y,pnt2[0].y+C.x, -pnt2[0].x+C.y);
          g.drawLine(C.x, C.y,pnt2[ndeg-1].y+C.x, -pnt2[ndeg-1].x+C.y);
        }
      for(i=0; i<= ndeg; i++) { x1[i] = Pnt1[i].x; y1[i] = Pnt1[i].y; w1[i] = Pnt1[i].z;}
      for(i=0; i< ndeg; i++) { //  投影面での座標
      x2[i] = Pnt1[i+1].x/Pnt1[i+1].z-Pnt1[i].x/Pnt1[i].z; 
	  y2[i] = Pnt1[i+1].y/Pnt1[i+1].z-Pnt1[i].y/Pnt1[i].z;}

	  // ２端点との２次元での距離比較
	  double x11=x1[0]/w1[0]; double y11=y1[0]/w1[0];
	  double x1n=x1[ndeg]/w1[ndeg]; double y1n=y1[ndeg]/w1[ndeg];
	  double dismin2=x11*x11 + y11*y11; // |PP0|
	  r2 = x1n*x1n + y1n*y1n; // |PPn|
      if(r2<dismin2) {dismin2=r2; xint[0]=1.;
	    xint[1]=x1n;xint[2]=y1n;}
	  else { xint[0]=0.;xint[1]=x11;xint[2]=y11;}

	  if(dismin2>dRmax*dRmax) { dismin2=dRmax*dRmax;xint[0]=-1.;}
      xint[3]=dismin2;
      nsplit = 0;   nint = 0;  n_clip = 0; nbar=0; 
      t1[0] = 0.0;  t1[1] = 1.0; 

      nint = d_ovrlp(ndeg,x1,y1,w1,t1,x2,y2,xint,g); // 前処理：距離でクリップ
	  if(nint<0) {
		  u_dist.u=xint[0];u_dist.v=0;
		g.setColor(Color.blue);
	    g.drawString(" no closest point within "+Rmax,20,44+12*Cnum);
	    return new Point(C.x,C.y); }

      nint = s_ovrlp(ndeg,x1,y1,w1,t1,x2,y2,xint,nint,g);

     g.setColor(Color.white);
	 if(xint[0]<0) {  //select_u=-1.; err のためサンプリング法を使用  
	    PntQ = clospnt2(ndeg,pnt1,xint1,g); 
		u_dist.u=xint1[0];u_dist.v=Math.sqrt(xint1[3]);
		g.setColor(Color.white);
	    g.drawString(" no closest point within "+Rmax,20,44+12*Cnum);
	    return new Point(C.x,C.y); }
	 else { 
	     double dist = Math.sqrt(xint[3]);
		 Point3 Pwrd = Bez[Cnum].eval(ndeg+1,xint[0]);
         u_dist.u =xint[0]; u_dist.v=dist;
		if(dist<DISmin) {DISmin=dist; select_u=xint[0];
           closestCurv=Cnum;
/*
		 if(Ncurve>0)
         g.drawString("Curve="+Cnum+" u= "+xint[0] +" d_min="+dist+
			" x="+Pwrd.y+" y="+Pwrd.z+" z="+Pwrd.x,20,44+12*Cnum);
		 else 
         g.drawString(" u= "+xint[0] +" d_min="+dist+
			" x="+Pwrd.y+" y="+Pwrd.z+" z="+Pwrd.x,20,44+12*Cnum);
*/
//		   closestPnt = new Point((int)xint[1]+C.x,(int)xint[2]+C.y);
		   if(xint[0]<0.01 || xint[0]>0.99) {
		   showStatus("u="+xint[0]+" distance="+Math.sqrt(xint[3])+
		  " Number of split="+ nsplit+" clip="+n_clip+" iterations="+nbar);}
		}
//		  " Number of split="+ nsplit+" clip="+n_clip+" u="+u_min);  
	 }
	 return new Point((int)xint[1]+C.x,(int)xint[2]+C.y); // 最近点の２次元座標を返す
    }
/* ***********サンプリイング法で距離を算出 **************: */
	public Point clospnt2(int ndeg,Point pnt1[],double xint[],Graphics g)   {
		Point  PntQ= new Point(0,0);
		double r2,u_min=0.;  int i;

		double dismin2=100000.;
		for( i=0; i <= ndiv; i++){
				double u = (double)i /(double)ndiv;
			    Point current = castljau(u,ndeg,pnt1);
				r2= (double)current.x*current.x+ current.y*current.y;
				if(r2<dismin2) {dismin2=r2; u_min=u;
					PntQ = new Point(current.x, current.y);}
        }
		xint[0]=u_min; 
		xint[1]=(double)PntQ.x; xint[2]=(double)PntQ.y;
		xint[3]=dismin2;
        return new Point(PntQ.x,PntQ.y);
    }
/* ***** 半径distの円の内接ベルトで曲線をクリップ **************: */
int d_ovrlp(int ndeg,double x1[],double y1[],double w1[],double t1[],double x2[],double y2[],
	double xint[],Graphics g) {
     double[] tlr = new double[2];
     double[] cf = new double[10];
	 double a,b,dist,span;
     int  i,iover,inega,iposi;

     tlr[0]=0.; tlr[1]=1.;
     span = 1./(double)ndeg;

     dist = Math.sqrt(xint[3]); 
     a = x1[ndeg]/w1[ndeg]-x1[0]/w1[0]; b= y1[ndeg]/w1[ndeg]-x1[0]/w1[0];
	 double dd =  Math.sqrt(a*a+b*b);
	 a =-a/dd; b= -b/dd; // 2端点ﾍﾞｸﾄﾙに一致する法線
     inega = 0; iposi = 0;
	 for(i=0;i<=ndeg; i++){
		 cf[i] = (a*x1[i]+b*y1[i])/w1[i] + dist;
              if(cf[i] > 0.0001) iposi = 1;
              else if(cf[i] < -0.0001) inega = 1;
     }
       //System.err.println(" cf="+cf[0]+", "+cf[1]+", "+cf[2]+", "+cf[3]);
       if(iposi == 0) { return (-1);} /* all negative means : dont't overlap */
       if(inega != 0) {             /*      (+ + + +) */
		   if(root3(ndeg, span, cf,tlr)<0) {return (-1);} 
            //System.err.println(" tlr="+tlr[0]+", "+tlr[1]);
         } 
     inega = 0; iposi = 0;
	 for(i=0;i<=ndeg; i++){
		 cf[i] = (-a*x1[i] - b*y1[i])/w1[i] + dist;
              if(cf[i] > 0.0001) iposi = 1;
              else if(cf[i] < -0.0001) inega = 1;
     }
     //System.err.println(" cf="+cf[0]+", "+cf[1]+", "+cf[2]+", "+cf[3]);
       if(iposi == 0) { return (-1);} /* all negative means : dont't overlap */
       if(inega != 0) {             /*      (+ + + +) */
		   if(root3(ndeg, span, cf,tlr)<0) {return (-1);} 
   //      System.err.println(" tlr="+tlr[0]+", "+tlr[1]);
         } 
    double tmp= t1[0];
    t1[0]=tmp+(t1[1]-tmp)*tlr[0]; t1[1]=tmp+(t1[1]-tmp)*tlr[1];
//    lr_split(ndeg,x1,y1,tlr[0],1.-tlr[1]); n_clip++ ;
    lr_Split(ndeg,x1,y1,w1,tlr[0],1.-tlr[1]); n_clip++ ;

    for(i=0; i < ndeg; i++) {
            x2[i] = x1[i+1]/w1[i+1]-x1[i]/w1[i]; y2[i] = y1[i+1]/w1[i+1]-y1[i]/w1[i];}
            //x2[i] = x1[i+1]-x1[i]; y2[i] = y1[i+1]-y1[i];}
	if(disp_con==0){
 	  g.setColor(Color.yellow);
	  for(i=0;i<ndeg;i++) g.drawLine((int)x1[i]+PntC.x,(int)y1[i]+PntC.y,
	  (int)x1[i+1]+PntC.x,(int)y1[i+1]+PntC.y);
    }
    return 0;
 }
/* *** サンプリング点での距離の最小値算出 **********************: */
public int caldismin(int ndeg,double x1[],double y1[],double w1[],
					 double xint[], int ndiv1,Graphics g){
     int nint=0;
     Point2 PntQ; 
	 double u,r2,u_min= 0.;

        double dismin2=xint[3];
		for(int i=1; i <= ndiv1; i++){
				 u = (double)i /(double)ndiv1;
//			    Point2 current = castljau2(u,ndeg,x1,y1);
		    Point3 current = castljau3(u,ndeg,x1,y1,w1);
				r2= (double)(current.x*current.x+ current.y*current.y)/current.z/current.z;
				if(r2<dismin2) {dismin2=r2; u_min=u; nint=1;
					xint[1]=current.x/current.z; xint[2]=current.y/current.z;}
        }
     if (nint>0) {xint[3] = dismin2; xint[0]=u_min;}
     return nint;
}
/* ********************************************************************: */
int s_ovrlp(int ndeg,double x1[],double y1[],double w1[],double t1[],double x2[],double y2[],
	double xint[],int nint,Graphics g) {
	 double[] x1l= new double[20]; double[] y1l= new double[20];
     double[] t1l= new double[2];
     double[] x1r= new double[20]; double[] y1r= new double[20];
	 double[] t1r= new double[2];
     double[] w1l= new double[20]; double[] w1r= new double[20];

     int  i,iover;

    iover = ovrlp(ndeg,x1,y1,w1,t1,x2,y2, xint,g);
	if(disp_con==0){// debug時 プロット
	  g.setColor(Color.green);
	  if(iover==0) g.setColor(Color.darkGray);
 	  if(iover==1) g.setColor(Color.white);
 	  for(i=1; i <ndeg; i++) {
       g.drawLine((int)y2[i-1]+PntC.x, -(int)x2[i-1]+PntC.y, 
			(int)y2[i]+PntC.x, -(int)x2[i]+PntC.y);}
       g.setColor(Color.darkGray);
       g.drawLine(PntC.x, PntC.y,(int)y2[0]+PntC.x, -(int)x2[0]+PntC.y);
       g.drawLine(PntC.x, PntC.y,(int)y2[ndeg-1]+PntC.x,-(int)x2[ndeg-1]+PntC.y);
    }
	if(iover == 1) { nint= nint +1; } // 解有り
    else if(iover == -1) { //曲線を２分割
          nsplit ++; if(nsplit>20) return 0;
//          subdiv(ndeg,x1,y1, x1l,y1l, x1r,y1r,0.5,t1,t1l,t1r);
//System.err.println(" ovlp;t0="+t1[0]+" t1="+t1[1]+"ndeg="+ndeg);
          subDiv(ndeg,x1,y1,w1,x1l,y1l,w1l,x1r,y1r,w1r,0.5,t1,t1l,t1r);
            for(i=0; i < ndeg; i++) {
            x2[i] = x1l[i+1]/w1l[i+1]-x1l[i]/w1l[i]; 
			y2[i] = y1l[i+1]/w1l[i+1]-y1l[i]/w1l[i];}
          nint=s_ovrlp(ndeg,x1l,y1l,w1l,t1l,x2,y2, xint,nint,g);
            for(i=0; i < ndeg; i++) {
            x2[i] = x1r[i+1]/w1r[i+1]-x1r[i]/w1r[i]; 
			y2[i] = y1r[i+1]/w1r[i+1]-y1r[i]/w1r[i];}
           nint=s_ovrlp(ndeg,x1r,y1r,w1r,t1r,x2,y2,xint,nint,g);
	}
    else if(iover == 2) { //再度よぶ
//		System.err.println(" iover=2");
          nint=s_ovrlp(ndeg,x1,y1,w1,t1,x2,y2, xint,nint,g);
	}
    return nint;
 }
/* ************************************************************: */
void plotrange(double t1[], int index,Graphics g)
{ int length =300,ybar=295,ybar1;

//    if(index==0) {nbar=0; return;}
	nbar++;
    if(nbar>30)	 System.err.println("Test:nbar="+nbar);
	if(idebug==0 && disp_mod==0) return;
	  ybar1=ybar-nbar;
	   g.setColor(Color.yellow);
	   if(nbar%10==0) g.setColor(Color.orange);
       g.drawLine(10,ybar1,10+length,ybar1);
 	   g.setColor(Color.red);
       g.drawLine(10+(int)(length*t1[0]),ybar1,10+(int)(length*t1[1]),ybar1);
}
/* ************************************************************: */
int ovrlp(int ndeg,double x1[],double y1[],double w1[],double t1[],
		  double x2[],double y2[],double xint[], Graphics g)
/*  Checks the overlap of the two curve segments
    return = 0  for no overlap (=no intersections)
          = -1 for inadequate clip (split each curve and try again)
          = 1  for intersection found                                     */
{
      double tmp1, tmp2, ur22;
      double  ur11, tl2,ur2,tl1, ur1;
	  double d1,d2,d,dd,dmin,dmax,a=0.,b=1.;
      int i,icount,jover=0;
	  double[] t2 = new double[2];
     
	 plotrange(t1,1,g); 
	  // test
	if(itest==0) {
	   int indQ = bmul1(x1,y1,w1,x2,y2,ndeg);
	   if(indQ==1) return 0;} // 極値はない
    if(itest==1) {
	   int indQ = bmul3(x1,y1,w1,x2,y2,ndeg,t2,t1,0,g);
	   if(indQ==1) return 0; // 極値はない
	   if(idebug==1) System.err.println("Q ="+indQ+" t="+t1[0]+", "+t1[1]
		   +" clip="+n_clip+" split="+nsplit);
        ur22 = 1. - t2[1];		  tl2=t2[0]; ur2=t2[1];
        tmp2 = t1[0] + tl2*(t1[1]-t1[0]);  tmp1 = t1[1] + ur22*(t1[0]-t1[1]);
           lr_Split(ndeg,x1,y1,w1,tl2,ur22); n_clip++ ;
		   t1[0]=tmp2; t1[1]=tmp1;
          if(idebug==1)
		  System.err.println(" aft g;lr_S; t10=" + t1[0]+",t1="+t1[1]+
			  " delt="+(tmp1-tmp2)+" tl2="+tl2+" ur2="+ur2);
		if((tmp1-tmp2)< TOLEL) {
		  if(disp_mod==0|| disp_mod==1) {g.setColor(Color.blue);
            if(tl2>0.02 || ur22>0.02) {
			    lr_Split(ndeg,x1,y1,w1,tl2,ur22); n_clip++ ;//clipping	
		    }
			 double r2 = (x1[0]*x1[0]+y1[0]*y1[0])/w1[0]/w1[0];
		     if(r2<xint[3]) {
		      xint[0]= t1[0]; xint[1]=x1[0]/w1[0];xint[2]=y1[0]/w1[0];
		      xint[3]= r2;
		     }
			 if(idebug==1)
			 System.err.println("find at g_hull;t1=" + t1[0]+" dt="+(tmp1-tmp2)
				 +" r2="+r2+" x="+xint[1]);
//			 if(disp_mod==0)
 //           g.drawLine(PntC.x, PntC.y,(int)xint[1]+PntC.x,(int)xint[2]+PntC.y);
		   return 1; // solved 
        }
        if(t1[1]-t1[0]>0.3) {
          // System.err.println(" t1-t0="+(t1[1]-t1[0]));
             for(int k=0;k<ndeg;k++)
             g.drawLine(PntC.x+(int)x1[k],PntC.y+(int)y1[k], 
			 PntC.x+(int)x1[k+1],PntC.y+ (int)y1[k+1]); }
			return -1;}
	       for(i=0; i < ndeg; i++) {
            x2[i] = x1[i+1]/w1[i+1]-x1[i]/w1[i]; 
			y2[i] = y1[i+1]/w1[i+1]-y1[i]/w1[i];}
		   if(tl2>0.1 || ur22>0.1) {
//			 System.err.println("return 2 g_hull;t1=" +t1[0]+"t2="+t1[1]+"tl2="+tl2);
		    return 2;}
	 }

	  d1 = Math.sqrt(x2[0]*x2[0]+y2[0]*y2[0]); // 端点の長さ
      a=x2[0]/d1;  b= y2[0]/d1;

      double d0 = (x1[0]*a + y1[0]*b)/w1[0];  
  	  d2 = Math.sqrt(x2[ndeg-1]*x2[ndeg-1]+y2[ndeg-1]*y2[ndeg-1]);
      double aa=x2[ndeg-1]/d2; double bb= y2[ndeg-1]/d2;
      double dn= (x1[ndeg]*aa + y1[ndeg]*bb)/w1[ndeg];
	  dd = (aa*b - bb*a);

	  int inv=1;
      if(dd<0) { dd = -dd; a=-a; b=-b;inv=-1;}
	  dmin =0.; dmax=dd;
	  for(i=1;i<ndeg-1;i++) {
       d2 = Math.sqrt(x2[i]*x2[i]+y2[i]*y2[i]);
       d= (x2[i]*b - y2[i]*a)/d2;
	   if(d<dmin) dmin = d;  if(d>dmax) {dmax = d;
	   if(inv>0) inv=i;
	   else inv=-i;}
      }
//      if(dmax>0.65 || dmin<-0.65 ) { 
      if(dmax>0.85 || dmin<-0.85 ) { 
		  if(t1[1]-t1[0]>0.1) return -1;} //分割が必要
	if(itest==0) {
	  int iposi=0; int inega=0;
      for(i=0;i<=ndeg;i++) {
       d= -(x1[i]*a + y1[i]*b)/w1[i];
	   if(d>0.0001) iposi=1;
	   if(d<0.0001) inega=1;
      }
      if(iposi==0 || inega==0) {
  	    int iposi1=0; int inega1=0;
        for(i=0;i<=ndeg;i++) {
         d= -(x1[i]*aa + y1[i]*bb)/w1[i];
	     if(d>0.0001) iposi1=1;
	     if(d<0.0001) inega1=1;
        }
       if(iposi1==0 || inega1==0) { 
//		   System.err.println(" kainashi? ");
		    return 0;  } //解無し？
	  }
  }
//      if(dmax>0.30 || dmin<-0.30) {  return -1;} //分割必要
      if(dmax>0.70 || dmin<-0.70) { 
		  if(t1[1]-t1[0]>0.1) return -1;} //分割必要
 
	  if(t1[1]-t1[0]<0.035) { // 無理やりサンプリング法で判定
	   int ndiv1 = (int)((t1[1]-t1[0])/0.005)+1; if(ndiv1>20) ndiv1=20;
	   int nint= caldismin(ndeg,x1,y1,w1,xint,ndiv1,g);
       if(nint>0) { xint[0]= (t1[1]-t1[0])*xint[0]+t1[0];
	      System.err.println(" caldismin u="+xint[0]+" x="+xint[1]+" r2="+xint[3]+
			  " ndiv="+ndiv1+" t0="+t1[0]+" t1="+t1[1]);
	   	  if(disp_mod==0){
			g.setColor(Color.yellow);
          g.drawLine(PntC.x, PntC.y,(int)xint[1]+PntC.x,(int)xint[2]+PntC.y);}
	      return 1;
	   }
	  }
      if(t1[1]-t1[0]<0.0001) {System.err.println(" t1<0.0005??");return 0;}
      t2[0]=0.; t2[1]=1.; iproces=0;
      for(icount=0; icount< 20; icount ++) {
          iproces++;         if(iproces>1)return 0;
/*  curve-x1 is clip by curve-x2 */
		jover = hull(ndeg,x1,y1,w1,x2,y2,a,b,aa,bb,t2,g); //ウエッジで分割範囲を判定
	    if(jover < 0) return -1;
        if(jover == 1) return -1;
        if(jover == 2) return -1;
        ur22 = 1. - t2[1];		  tl2=t2[0]; ur2=t2[1];
        tmp2 = t1[0] + tl2*(t1[1]-t1[0]);  tmp1 = t1[1] + ur22*(t1[0]-t1[1]);
        t1[1] = tmp1;  t1[0] = tmp2;
		if((tmp1-tmp2)< TOLEL) {
		     if(disp_mod==0) {g.setColor(Color.white);
             for(int k=0;k<ndeg;k++)
             g.drawLine(PntC.x+(int)x1[k],PntC.y+(int)y1[k], 
			 PntC.x+(int)x1[k+1],PntC.y+ (int)y1[k+1]); }

			 lr_Split(ndeg,x1,y1,w1,tl2,ur22); n_clip++ ;//clipping	
			 double r2 = (x1[0]*x1[0]+y1[0]*y1[0])/w1[0]/w1[0];
		   if(r2<xint[3]) {
		    xint[0]= t1[0]; xint[1]=x1[0]/w1[0];xint[2]=y1[0]/w1[0];
		    xint[3]= r2;
      	    if(disp_mod==0) {g.setColor(Color.green);
            g.drawLine(PntC.x, PntC.y,(int)xint[1]+PntC.x,(int)xint[2]+PntC.y);}
		   }
	        if(disp_mod==0)
			 System.err.println("find at hull;t1=" + t1[0]+" dt="+(tmp1-tmp2)
				 +" r2="+r2+" x="+xint[1]+" tl2="+tl2+" ur2="+ur2);
		   return 1; // solved 
        }
       if(tl2> 0.02 || ur2 < 0.98) {
           lr_Split(ndeg,x1,y1,w1,tl2,ur22); n_clip++ ;
//		   showStatus(" aft hul1;Number of split="+ nsplit+" Number of clip="+n_clip);
//		  System.err.println(" aft lr_S; t10=" + t1[0]+",t1="+t1[1]+
//			  "delt="+(tmp1-tmp2)+" tl2="+tl2+" ur2="+ur2);
       if(t1[1]-t1[0]<0.05){
		   //  test 無理やりサンプリング法で判定
	     int ndiv1 = (int)((t1[1]-t1[0])/0.005)+1; if(ndiv1>25) ndiv1=25;
	     int nint= caldismin(ndeg,x1,y1,w1,xint,ndiv1,g);
         if(nint>0) { xint[0]= (t1[1]-t1[0])*xint[0]+t1[0];
//	       System.err.println(" aft lr_S;caldismin u="+xint[0]+" x="+xint[1]+" r2="+xint[3]+
//			  " ndiv="+ndiv1);
	   	   if(disp_mod==0){
		 	g.setColor(Color.yellow);
            g.drawLine(PntC.x, PntC.y,(int)xint[1]+PntC.x,(int)xint[2]+PntC.y);}
	        return 1;
	       }
	   }
	   // end of test
	   	  if(disp_mod==0){  g.setColor(Color.pink);
             for(int k=0;k<ndeg;k++)
             g.drawLine(PntC.x+(int)x1[k],PntC.y+(int)y1[k], 
			 PntC.x+(int)x1[k+1],PntC.y+ (int)y1[k+1]); } 
	      if(disp_mod==0)
		  System.err.println(" ; t10=" + t1[0]+",t1="+t1[1]+
			  "delt="+(tmp1-tmp2));
		  return -1; //再分割が必要
       }
	// これ以後は変更予定（未完成）
/*  curve-1 is clip by curve-2 */
     //  jover = hull(ndeg,x1,y1,g);
	   jover =0;
//	   System.err.println(" hull-2:jover="+jover);
       if(jover < 0) return 0;
       if(jover == 1) return (-2);
//	   System.err.println(" hull-2:tl1=" + tl1+",ur1="+ur1+" jover="+jover);
/*
       if((ur2 - tl2 > 0.80)) return (-1) ;

        ur11 = 1-ur1;
        tmp1 = t1[0] + tl1*(t1[1]-t1[0]); tmp2 = t1[1] + ur11*(t1[0]-t1[1]);
       if( (ur2 - tl2 > 0.98) && icount>0 ) {
           t1[1] = tmp2;  t1[0] = tmp1;
           lr_split(ndeg,x1,y1,tl1,ur11); n_clip++ ;
//		   System.err.println(" x1:split; tl1=" + tl1+",ur1="+ur1);
//		 showStatus(" Number of split="+ nsplit+" Number of clip="+n_clip);
		 if(disp_con==0) {
		   g.setColor(Color.gray);
           g.drawLine((int)x1[0], (int)y1[0], (int)x1[ndeg], (int)y1[ndeg]);
		 } 
           return (-3) ;      }

       if((tmp2-tmp1)<= TOLEL) {
           if((t1[1]-t1[0])<TOLEL ) {
             t1[1] = tmp2;  t1[0] = tmp1;
             return (1) ;
           }
       }
 //      if(tl1<0.01 && ur1>0.99)  return (-3) ;   
       t1[1] = tmp2;  t1[0] = tmp1;
       lr_split(ndeg,x1,y1,tl1,ur11); n_clip++; 
	   */
//	   showStatus(" Number of split="+ nsplit+" Number of clip="+n_clip);
      //   for(int k=0;k<ndeg;k++)
        // g.drawLine((int)x1[k], (int)y1[k], (int)x1[k+1], (int)y1[k+1]);} 
       }
      return 0;
  }
/* ******** ウエッジで曲線をクリップ *********************: */
int hull(int ndeg,double x[],double y[],double w[], double x2[],
		 double y2[],double a,double b,double aa,double bb,
		 double t2[],Graphics g) {
     int  i, il, iposi,k, inega,inega1=0;
     double x0,y0,x3,y3,tmp,tmp1,span;
     double[] cf = new double[10];
     double[] tlr = new double[2];
     double tl, ur;
     int nside = 2;

     span = 1./(double)ndeg;
     tl = 0.; ur = 1.; //TU= new Vec2(tl,ur);
     tlr[0]=0.; tlr[1]=1.;

     for(il=1; il <= nside; il++) {
            inega = 0; iposi = 0;
      if(il == 1) {
            for(k=0; k <= ndeg; k++) {
                cf[k] = -(x[k]*a+x[k]*b)/w[k];
              if(cf[k] > 0.0001) iposi = 1;
              else if(cf[k] < -0.0001) inega = 1;}
			 inega1=inega;
             }
      else if(il == 2) {
            for(k=0; k<= ndeg; k++) {
              cf[k] =  -(x[k]*aa+x[k]*bb)/w[k];
              if(cf[k] > 0.0001) iposi = 1;
              else if(cf[k] < -0.0001) inega = 1;
          }
			if(inega==0 && inega1==0) return 2;
         }
 //   System.err.println(" il="+il+" cf="+cf[0]+", "+cf[1]+", "+cf[2]+", "+cf[3]);
       if(iposi == 0 && inega==1) { return (1);} /* all negative means : dont't overlap */
       if(inega != 0 && iposi!=0) {     /*      (+ + + +) */
		   if(root3(ndeg, span, cf,tlr)<0) { return (-1);} 
   //        System.err.println(" tlr="+tlr[0]+", "+tlr[1]);
		   t2[0]=tlr[0];t2[1]=tlr[1];
         } 
     }   /* loop for il */
     return (0);
}
/* *********************** ２項係数の計算 **************************: */
public void binit(int maxdeg){
	double mum,bc,den;

	for(int i=0;i<maxdeg;i++){
		 mum =(double)i;
		bc = 1.;
		den =1.;
		for(int j=0; j<=i;j++){
			binom[i][j]= bc;
			bc = bc*mum/den;
			mum = mum -1.0; den = den +1.;
//	       System.err.println(i+" bin="+binom[i][j]);
		}
	}
}
/* ************************************************************: */
public int bmul2(double x1[],double y1[],double x2[],double y2[],int ndeg){

	  double[] x3 = new double[30];
	  double[] y3 = new double[30];

      int ndeg3 =bmul(x1,ndeg,x2,ndeg-1,x3);
       ndeg3 =bmul(y1,ndeg,y2,ndeg-1,y3);

	   int nega=0, posi=0;
	   for(int i =0;i<=ndeg3; i++) {
		  x3[i] += y3[i];
		  if(x3[i]>0) posi=1;
		  if(x3[i]<0) nega=1;
//	       System.err.println(i+" x3="+x3[i]);
	  }
	   if(nega==0 || posi==0) return 1;
	   return 0;
}
/* ************************************************************: */
public int bmul1(double x1[],double y1[],double w1[],double x2[],double y2[],int ndeg){
      int i,i1,i2,i3,n3;
	  double[] x3 = new double[20];

	  int n2 =ndeg-1;
 	  n3 = ndeg+n2;
	  x3[0]= (x1[0]*x2[0]+y1[0]*y2[0])/w1[0];
	  x3[n3]=(x1[ndeg]*x2[n2]+y1[ndeg]*y2[n2])/w1[ndeg];
	  if(x3[0]>0.001 && x3[n3]<-0.001) return 0;
	  if(x3[0]<-0.001 && x3[n3]>0.001) return 0;

	  for(i=1; i<n3; i++) x3[i]=0.;

	for(i1=0;i1<=ndeg;i1++) {
		for(i2=0; i2<ndeg; i2++) {
			i3 = i1+i2;
			if(i3>0 && i3<n3) 
			x3[i3] += 
		binom[ndeg][i1]*binom[n2][i2]*(x1[i1]*x2[i2]+y1[i1]*y2[i2])/w1[i1];
		}
	}

	   if(x3[0]>0.0) {
	   for( i =1;i<n3; i++) {
		  if(x3[i]<-0.001) return 0;
	   //    System.err.println(i+" x3="+x3[i]);
	  }
	   }
	   else{
		  for( i =1;i<n3; i++) {
		  if(x3[i]>0.001) return 0;
	    //   System.err.println(i+" x3="+x3[i]);
	      }
	   }
	   return 1;
}
/* ************************************************************: */
public int bmul3(double x1[],double y1[],double w1[],double x2[],
				 double y2[],int ndeg, double t2[],double t1[],
				 int mode, Graphics g){
      int i,i1,i2,i3,n3;
	  double[] x3 = new double[30];
	 // double[] t2 = new double[2];

	  int n2 =ndeg-1;
 	  n3 = ndeg+n2;
	  x3[0]= (x1[0]*x2[0]+y1[0]*y2[0])/w1[0];
	  x3[n3]=(x1[ndeg]*x2[n2]+y1[ndeg]*y2[n2])/w1[ndeg];

	  for(i=1; i<n3; i++) x3[i]=0.;
	for(i1=0;i1<=ndeg;i1++) {
		for(i2=0; i2<ndeg; i2++) {
			i3 = i1+i2;
			if(i3>0 && i3<n3) 
			x3[i3] += 
		binom[ndeg][i1]*binom[n2][i2]*(x1[i1]*x2[i2]+y1[i1]*y2[i2])/w1[i1];
		}
	}
//		System.err.println("Bef q_hull: x0="+x3[0]+" xn3="+x3[n3]);

	   if(x3[0]>0.0) {
		   int ine = 0;
	   for( i =1;i<=n3; i++) {
		   x3[i]=x3[i]/binom[n3][i];
		  if(x3[i]<-0.001) ine = 1;
//	       System.err.println(i+" x3="+x3[i]);
	  }
	   if(ine==0) return 1;
	   }
	   else{
		   int ipo = 0;
		  for( i =1;i<=n3; i++) {
		   x3[i]=x3[i]/binom[n3][i];
		  if(x3[i]>0.001) ipo=1;
//	       System.err.println(i+" x3="+x3[i]);
	      }
		  if(ipo==0) return 1;
	   }
if(idebug==1) {
	   int yaxis=170; double sc=0.01;double length=200.;
	   g.setColor(Color.white);
	   g.drawLine(10,yaxis,10+(int)length,yaxis);
	   nbar++;
	   g.setColor(Color.red);
	   if(t1[1]-t1[0]<0.3) {sc=0.1; g.setColor(Color.green);}
       if(mode>0) g.setColor(Color.white);

	   for( i=0; i < n3; i++){
		   double length1 =length*(t1[1]-t1[0]);
		   int offset = (int)(length*t1[0])+10;
		   int ix1 = offset+(int)(((double)i/(double)n3)*length1);
           int ix2 = offset+(int)(length1*(double)(i+1)/(double)n3);
		   g.drawLine(ix1,(int)(x3[i]*sc)+yaxis,ix2,(int)(x3[i+1]*sc)+yaxis);
		}
}
	   int jover = q_hull(n3,x3, t2);
	    if(jover < 0) return 1;
        if(jover == 1) return 1;
        if(jover == 2) return 1;
        double ur22 = 1. - t2[1];		  
		double tl2=t2[0]; double ur2=t2[1];
//		System.err.println("q_hull: tl2="+tl2+" ur2="+ur2);
       return 0;
}
/* ************************************************************: */
public int bmul(double p1[],int n1,double p2[], int n2,double p3[]){
int i,i3,n3;

	n3 = n1 +n2;
    for( i=0; i<=n3; i++) p3[i]=0.;

	for(int i1=0;i1<=n1;i1++) {
		for(int i2=0; i2<=n2; i2++) {
			i3 = i1+i2;
			p3[i3]= p3[i3]+binom[n1][i1]*binom[n2][i2]*p1[i1]*p2[i2];
		}
	}

	for(i=1;i<=n3; i++) {
		p3[i]= p3[i]/binom[n3][i];
//	       System.err.println(i+" p3="+p3[i]);
	}
	return n3;
}
int q_hull(int ndeg,double ck[], double t2[]) {
     int  i, il, iposi,k, inega,inega1=0;
     double x0,y0,x3,y3,tmp,tmp1,span;
     double[] cf = new double[20];
     double[] tlr = new double[2];
     double tl, ur;
     int nside = 2;

     span = 1./(double)ndeg;
     tl = 0.; ur = 1.; //TU= new Vec2(tl,ur);
     tlr[0]=0.; tlr[1]=1.;

     for(il=1; il <= nside; il++) {
            inega = 0; iposi = 0;
      if(il == 1) {
            for(k=0; k <= ndeg; k++) {
				cf[k] = ck[k];
              if(cf[k] > 0.0001) iposi = 1;
              else if(cf[k] < -0.0001) inega = 1;}
			 inega1=inega;
             }
      else if(il == 2) {
            for(k=0; k<= ndeg; k++) {
              cf[k] =  -ck[k];
              if(cf[k] > 0.0001) iposi = 1;
              else if(cf[k] < -0.0001) inega = 1;
          }
			if(inega==0 && inega1==0) return 2;
         }
//	  System.err.println(" q:il="+il+" cf="+cf[0]+", "+cf[1]+", "+cf[2]+", "+cf[3]
//		  +", "+cf[4]+", "+cf[5]);
       if(iposi == 0 && inega==1) { return (1);} /* all negative means : dont't overlap */
       if(inega != 0 && iposi!=0) {     /*      (+ + + +) */
		   if(root3(ndeg, span, cf,tlr)<0) { return (-1);} 
 //          System.err.println(" tlr="+tlr[0]+", "+tlr[1]);
		   t2[0]=tlr[0];t2[1]=tlr[1];
         } 
     }   /* loop for il */
     return (0);
}

/* ************ 解の存在範囲の抽出 *********************: */
public int root3(int nmdeg, double span, double cf[], double tlr[]) {
      double t, tmin, tmax;
      int i,j;

      if((cf[0] > 0.) && (cf[nmdeg] < 0.) ) /* (+ ? ? -) > right cut [0,tmax] */
      {     tmax = 0. ;  /* do 310 */
            for (i=0; i < nmdeg; i++) {
              if(cf[i] > 0.) {
                for(j= i+1 ; j <= nmdeg; j++) {
                  if(cf[j] < cf[i]) {
                    if(cf[j] < 0.) {
                    t = (cf[i]/(cf[i]-cf[j])*(double)(j-i)+i)*span;
                    if(t > tmax) tmax = t;
                    }
                  }
                }
              }
             }
            if(tmax != 0.) {
              if(tlr[0] > tmax+.000001) return (-1)  ;
              if(tlr[1] > tmax) tlr[1] = tmax;}
       } 
       else if( (cf[0] < 0.) && (cf[nmdeg] > 0.) ) 
/*                            (- ? ? +) > left cut  [tmin,1.] */
       {    tmin = 1. ;  /* do 320 */
            for(i=0; i < nmdeg; i++) {
             if(cf[i] < 0.) {
              for(j=i+1; j <= nmdeg ; j++) {
                 if(cf[j] > cf[i]) {
                     if(cf[j] > 0.) {
                     t = (cf[i]/(cf[i]-cf[j])*(double)(j-i)+i)*span ;
                     if(t < tmin) tmin = t ;
                     }
                 }
               }
             }
            }
            if(tmin != 1.) {
                  if(tlr[1] < tmin-0.000001) return (-1) ;
                  if(tlr[0] < tmin) tlr[0] = tmin;   }
        }
        else if( (cf[0]> 0) && (cf[nmdeg] > 0.) ) { }
       /*               split 1/2          */
        else if( (cf[0] < 0.) && (cf[nmdeg] < 0.)) 
/*                         (- ? ? -) > both side cut    [tmin,tmax] */
        {   tmax = 0. ;       /* do 330 */
            for(i=1; i < nmdeg; i++) {
                for(j=i+1; j <= nmdeg; j++) {
                    if( (cf[i] > 0.) && (cf[j] < 0.)) {
                    t = (cf[i]/(cf[i]-cf[j])*(double)(j-i)+i)*span ;
                    if(t > tmax) tmax = t ;
                    }
                }
             }
            if(tmax != 0) {
                   if(tlr[0] > tmax+.000001) return (-1)  ;
                   if(tlr[1] > tmax) tlr[1] = tmax ; }
           tmin = 1. ;       /* do 340 */
           for(i=0; i < nmdeg-1; i++)    {
              for(j=i+1; j < nmdeg; j++) {
                 if((cf[i] < 0) && (cf[j] > 0.)) {
                 t = (cf[i]/(cf[i]-cf[j])*(double)(j-i)+i)*span ;
                 if(t < tmin) tmin = t;     
               }
               }
            }
            if(tmin != 1.) {
               if(tlr[1] < tmin-0.000001) return (-1) ;
               if(tlr[0] < tmin) tlr[0] = tmin ;  }
       }
         return (0) ;
}
/* ***********サンプリング法で距離を算出 **************: */
	public Point clospnt1(int ndeg,Point pnt[], Point C, Graphics g)   {
		Point  PntQ= new Point(0,0); Point[] pnt1 = new Point[20];
		double r2,u_min=0.;  int i;
        double[] xint1= new double[4];

		for( i=0;  i <=ndeg; i++) {
		pnt1[i] =  new Point(pnt[i].x-C.x,pnt[i].y-C.y);}
	    PntQ = clospnt2(ndeg,pnt1,xint1,g);  
  //    showStatus(" min_distance="+Math.sqrt(dismin2)+" u="+u_min);
      return new Point(PntQ.x+C.x,PntQ.y+C.y);
    }
/* *************** 描画関数（塗りつぶし）*****************: */
	/*
    public void Bezscan(int ndeg,Point pnt[],Graphics g) {
	 int i,ix,nint=0;
	 double[] xint = new double[4]; // xint[0]=u,xint[1]=x,xint[2]=y,xint[3]=distance
     double[] xint1= new double[4];
	 Point C;
	 Point  PntQ= new Point(0,0);
	 Point[] pnt1 = new Point[11];	Point[] pnt2 = new Point[10];
	 double r2,u_min=0.; 
     double[] t1= new double[2];
     double[] x1= new double[11]; double[] y1= new double[11];
     double[] x2= new double[10]; double[] y2= new double[10];
	 double[] w1= new double[10];
   
	 if(Rmax>40) Rmax=40;
	 int ymax=0; int ymin=320;
	 int xmax=0; int xmin=320;
	 for(i=0;i<=ndeg;i++){
		 if(pnt[i].y < ymin) ymin=pnt[i].y;
		 if(pnt[i].y > ymax) ymax=pnt[i].y;
		 if(pnt[i].x < xmin) xmin=pnt[i].x;
		 if(pnt[i].x > xmax) xmax=pnt[i].x;
	}
	 System.err.println(" ymin="+ymin+", ymax="+ymax);
	 System.err.println(" xmin="+xmin+", xmax="+xmax);
     int ymin1 = ymin-45; if(ymin1<20) ymin1=20;
     int xmin1 = xmin-55; if(xmin1<1) xmin1=1;
     int ymax1 = ymax+60; if(ymax1>349) ymax1=349;
     int xmax1 = xmax+55; if(xmax1>338) xmax1=338;

	 for(int iy=ymin1;iy<ymax1;iy=iy+3) {
		 for( ix=xmin1;ix<xmax1;ix=ix+3){
          C = new Point(ix,iy);
		  for( i=0;  i <=ndeg; i++) { // 平行移動
		  pnt1[i] =  new Point(pnt[i].x-C.x,pnt[i].y-C.y);}
		  for( i=0;  i <ndeg; i++) { //ホドグラフ計算
		  pnt2[i] =  new Point(pnt[i+1].x-pnt[i].x,pnt[i+1].y-pnt[i].y); }
          PntC = new Point(C.x,C.y);
          for(i=0; i<= ndeg; i++) { 
			  x1[i] = pnt1[i].x; y1[i] = pnt1[i].y;
			  w1[i]=1.;
		  }
          for(i=0; i< ndeg; i++) {
          x2[i] = pnt1[i+1].x-pnt1[i].x; y2[i] = pnt1[i+1].y-pnt1[i].y;}
	    // ２端点との距離比較
	      double dismin2=x1[0]*x1[0] + y1[0]*y1[0];
	      r2 = x1[ndeg]*x1[ndeg] + y1[ndeg]*y1[ndeg];
          if(r2<dismin2) {dismin2=r2; xint[0]=1.;
	       xint[1]=x1[ndeg];xint[2]=y1[ndeg];}
	       else { xint[0]=0.;xint[1]=x1[0];xint[2]=y1[0];}
     
	      if(dismin2>Rmax*Rmax) 
		       { dismin2=Rmax*Rmax;xint[0]=-1.;}
           xint[3]=dismin2;
           nsplit = 0;   nint = 0;  n_clip = 0;nbar=0;  
           t1[0] = 0.0;  t1[1] = 1.0; 
           nint = d_ovrlp(ndeg,x1,y1,w1,t1,x2,y2,xint,g); // 前処理：距離でクリップ
		   if(nint<0) {break;}

           nint = s_ovrlp(ndeg,x1,y1,w1,t1,x2,y2,xint,nint,g);
           //  g.setColor(Color.white);
	      if(xint[0]<0) {  // select_u=-1.;
	       PntQ = clospnt2(ndeg,pnt1,xint1,g);  
            g.setColor(Color.gray);
           if(xint1[3]<Rmax*Rmax) g.setColor(Color.green);
           if(xint1[3]<Rmax*Rmax/4) g.setColor(Color.pink);
  		  	g.fillOval(C.x-2, C.y-2, 4, 4);
	      }
	      else { select_u=xint[0];
            g.setColor(Color.yellow);
           if(xint[3]<Rmax*Rmax/4) g.setColor(Color.pink);
           if(xint[3]<Rmax*Rmax/9) g.setColor(Color.red);
  		  	g.fillOval(C.x-2, C.y-2, 4, 4);
	      }
		 }
	 }
	 }
	 */
/* ********* 曲線を２分割 *******************************: */
public void subdiv(int ndeg,double x[],double y[], 
			double xl[],double yl[], double xr[],
			double yr[],double t,double t0[],double tl[],double tr[]) {
     double[][] bx= new double[20][20];
     double[][] by= new double[20][20];
     int i,k, ii;

       for(i=0;i <= ndeg; i++) {
          bx[0][i] = x[i];  by[0][i] = y[i];
        }
       for(k=0; k <= ndeg; k++) {
         for(i=1; i <= ndeg; i++) {
             bx[k+1][i] = bx[k][i-1] + t*(bx[k][i] -bx[k][i-1]);
             by[k+1][i] = by[k][i-1] + t*(by[k][i] -by[k][i-1]);
          }
       }

      for(i=0;i <= ndeg; i++) {
         xl[i] = bx[i][i]; yl[i] = by[i][i]; 
         ii = ndeg-i;
         xr[i] = bx[ii][ndeg]; yr[i] = by[ii][ndeg]; 
      }
        tl[0] = t0[0]; tl[1] = t0[0]+(t0[1]-t0[0])*t;
        tr[0] = tl[1]; tr[1] = t0[1];
}
/* ********* 曲線を２分割 *******************************: */
public void subDiv(int ndeg,double x[],double y[], double w[],
			double xl[],double yl[], double wl[],double xr[],
			double yr[],double wr[],double t,double t0[],double tl[],double tr[]) {
     double[][] bx= new double[10][10];
     double[][] by= new double[10][10];
     double[][] bw= new double[10][10];

     int i,k, ii;

       for(i=0;i <= ndeg; i++) {
          bx[0][i] = x[i];  by[0][i] = y[i];  bw[0][i] = w[i];
        }
       for(k=0; k <= ndeg; k++) {
         for(i=1; i <= ndeg; i++) {
             bx[k+1][i] = bx[k][i-1] + t*(bx[k][i] -bx[k][i-1]);
             by[k+1][i] = by[k][i-1] + t*(by[k][i] -by[k][i-1]);
             bw[k+1][i] = bw[k][i-1] + t*(bw[k][i] -bw[k][i-1]);
          }
       }

      for(i=0;i <= ndeg; i++) {
         xl[i] = bx[i][i]; yl[i] = by[i][i]; wl[i] = bw[i][i];
         ii = ndeg-i;
         xr[i] = bx[ii][ndeg]; yr[i] = by[ii][ndeg]; wr[i] = bw[ii][ndeg];
      }
        tl[0] = t0[0]; tl[1] = t0[0]+(t0[1]-t0[0])*t;
        tr[0] = tl[1]; tr[1] = t0[1];
}
/* ********** 曲線描画 ***********************************: */
/*	public void DrawCurve(int ndeg,Point pnt[], Graphics g)   {
		Point old = pnt[0];
		for(int i=1; i <= ndiv; i++){
				double u = (double)i /(double)ndiv;
			    Point current = castljau(u,ndeg,pnt);
				g.drawLine(old.x,old.y,current.x,current.y);
				old = current;	}
    }
	*/
/* ************* 曲線上の座標評価 (integer)****************************: */
	public Point castljau(double u, int ndeg, Point pnt[])	{
       double[] wx = new double[10]; double[] wy = new double[10];

       for (int i=0; i <=ndeg; i++){
        wx[i] = (double)pnt[i].x;  wy[i] = (double)pnt[i].y;}
       for (int m=1; m <= ndeg; m++) {
        for (int j=0; j <= ndeg - m; j++) {
          wx[j] +=(wx[j+1]-wx[j])*u ;
          wy[j] +=(wy[j+1]-wy[j])*u;  }
       }
       return new Point((int)wx[0],(int)wy[0]);    
    }
/* ***************曲線上の座標評価 （Real） ************: */
	public Point2 castljau2(double u, int ndeg, double x1[], double y1[]){
       double[] wx = new double[10];  double[] wy = new double[10];

       for (int i=0; i <=ndeg; i++){
        wx[i] = x1[i];  wy[i] = y1[i];}
       for (int m=1; m <= ndeg; m++) {
        for (int j=0; j <= ndeg - m; j++) {
          wx[j] +=(wx[j+1]-wx[j])*u ;
          wy[j] +=(wy[j+1]-wy[j])*u;  }
       }
       return new Point2(wx[0],wy[0]);    
    }
/* ***************曲線上の座標評価 （Real） ************: */
	public Point3 castljau3(double u, int ndeg, double x1[], 
		double y1[], double w1[]){
       double[] wx = new double[10];  
	   double[] wy = new double[10];
	   double[] ww = new double[10];

       for (int i=0; i <=ndeg; i++){
        wx[i] = x1[i];  wy[i] = y1[i]; ww[i] = w1[i];}
       for (int m=1; m <= ndeg; m++) {
        for (int j=0; j <= ndeg - m; j++) {
          wx[j] +=(wx[j+1]-wx[j])*u ;
          wy[j] +=(wy[j+1]-wy[j])*u; 
          ww[j] +=(ww[j+1]-ww[j])*u; 
		}
       }
       return new Point3(wx[0],wy[0],ww[0]);    
    }
/* ******************** 曲線の左クリップ *********************************: */
	public void lsplit(int ndeg,double x[],double y[],double ur){
     int i,j;

      for(i=1; i<= ndeg; i++)
         for(j=ndeg; j>=i; j--) {
         x[j] = x[j] + (x[j-1] - x[j])*ur;
         y[j] = y[j] + (y[j-1] - y[j])*ur;
         }
 }
/* ***************** 曲線の右クリップ ****************************: */
 void rsplit(int ndeg,double x[],double y[],double tl)  {
     int i,j;

     for(i=1; i <= ndeg; i++)  {
         for(j=0; j <= ndeg-i; j++) {
         x[j] = x[j] + (x[j+1] - x[j])*tl;
         y[j] = y[j] + (y[j+1] - y[j])*tl;
         }
     }
  }
/* *****************曲線の左右クリップ*************************: */
	public void lr_split(int ndeg,double x[],double y[],double tl,double ur){
     double tmp;

      rsplit(ndeg,x,y,tl);  if(ur == 0.) return;
      tmp = ur/(1.0-tl); lsplit(ndeg,x,y,tmp);
  }
/* *****************曲線の左右クリップ*************************: */
	public void lr_Split(int ndeg,double x[],double y[],double w[],
		double tl,double ur){
     double tmp;

      rSplit(ndeg,x,y,w,tl);  if(ur < 0.00005) return;
      tmp = ur/(1.0-tl); lSplit(ndeg,x,y,w,tmp);
  }
/* ******************** 曲線の左クリップ ****************************************: */
	public void lSplit(int ndeg,double x[],double y[],double w[],double ur){
     int i,j;

      for(i=1; i<= ndeg; i++)
         for(j=ndeg; j>=i; j--) {
         x[j] = x[j] + (x[j-1] - x[j])*ur;
         y[j] = y[j] + (y[j-1] - y[j])*ur;
         w[j] = w[j] + (w[j-1] - w[j])*ur;
         }
 }
/* ***************** 曲線の右クリップ ****************************: */
 void rSplit(int ndeg,double x[],double y[],double w[],double tl)  {
     int i,j;

     for(i=1; i <= ndeg; i++)  {
         for(j=0; j <= ndeg-i; j++) {
         x[j] = x[j] + (x[j+1] - x[j])*tl;
         y[j] = y[j] + (y[j+1] - y[j])*tl;
         w[j] = w[j] + (w[j+1] - w[j])*tl;
         }
     }
  }
/* **************************************************************: */
    public boolean mouseDown(Event evt, int x, int y) {
		if(n < Nc && addcirc==0){ 
			   p[n] = new Point(x,y);	n++;
		}else{	n=0; }

        if(addcirc==1 && idrag==0) { 
		if(y<260) {ncirc++; inc_circ=1;
	        metaball[ncirc] = new Metaball((double)x,(double)y,0.,0.,1.);
	    	showStatus(" set Point = "+ncirc+" x="+x+" y="+y);
		}
		}

		prev = new Point(x, y);
		repaint();
		return true;
    }
/* ---------------------------------------------------------------------- */
	public boolean mouseDrag(Event ev, int x, int y) {
		 int ix,iy,id;

        if(findp==1) {
         double[] x1 = new double[10];double[] y1=new double[10];

		 ix = x-prev.x;	iy = y-prev.y;
		 if(ix==0 && iy==0) return true;
		 if( iy>295) return true;

		 for(int i=0;i<Nc-1;i++) {x1[i]=0.;y1[i]=0.;}
         if(select_u<0.) return true;
		  if(select_u<1./(double)Nc*0.2) id=0;
		  if(select_u>1-1./(double)Nc*0.2) id=Nc-1;
          id = (int)(select_u*(double)(Nc-1)+0.4);
		  if(id>Nc) id = Nc-1;
		  if(select_u<1./(double)Nc*0.2) id=0;
		  if(select_u>1-1./(double)Nc*0.2) id=Nc-1;
		  int iix = ix; int iiy=iy;	 x1[id]=1.;
		   Point2 current = castljau2(select_u,Nc-1,x1,y1);
		  if(current.x<0.25) current.x=0.25;
		  double xx = ((double)ix/current.x);
		  double yy = ((double)iy/current.x);

		  Point3 Dp = new Point3(0,xx,yy); 

//  System.err.println("Dp x="+Dp.x+" y="+Dp.y+" z="+Dp.z+"xx="+xx+"yy="+yy);
         //  xx = Dp.y; yy =Dp.z;
           Bez[closestCurv].P3[id] = new Point3(Bez[closestCurv].P3[id].x+Dp.x,
		   Bez[closestCurv].P3[id].y+Dp.y,Bez[closestCurv].P3[id].z+Dp.z);
		   Bez[closestCurv].p3[id] = tranfrm(Bez[closestCurv].P3[id]);
		   showStatus("Curv #="+closestCurv+" Moving Point= "+id+" u="+select_u+" scal="+(1./current.x));
           idrag=1;
		   repaint();
		   prev.move(x, y);
		}
		if(addcirc==1) {
		  ix = Math.abs(x-prev.x);
		  iy = Math.abs(y-prev.y);
		  id= ix; if(ix<iy) id=iy; 
          metaball[ncirc].r=(double)id;
          metaball[ncirc].r2=metaball[ncirc].r*metaball[ncirc].r;
    	  showStatus("Moving Point = "+ncirc+" r="+id);
          idrag=1;
		  repaint();
		}
		return true;
	}

/* **************************************************************: */
	public boolean mouseUp(Event ev, int x, int y) {
		prev = null;
        idrag=0;
		return true;
	}	
/* **************************************************************: */
	public boolean action(Event evt, Object obj) {
		if (evt.target == b1) { //new curve
			n=0;curve1=0;
			Ncurve=0;
			Rmax = Integer.parseInt(t1.getText());
			if(Nmesh>0) { mesh_move=0;findp=0;}
			repaint();
		} 
		if (evt.target == b2) { // add curve
			n=0;curve1=0;
//			System.err.println("New curve Nc="+Nc+" p1="+p1[0].x);
			Ncurve++;
			System.err.println(" Ncurve="+Ncurve);
			Rmax = Integer.parseInt(t1.getText());
			repaint();
		} 
		else if (evt.target == b3) { // debug button
			idebug= 1-idebug;
			repaint();
		} 
		else if (evt.target == b4) {
			if(disp_anime>0) { 
//			  System.err.println(" disp_anime="+disp_anime);
			  crc_anime() ;
			}
			else {
				if(disp_anime==0) disp_anime=-1;
				if(disp_anime==-1) disp_anime=1;
			for(int k=0; k<=Ncurve;k++) {
			for(int id=0;id<Nc;id++) {
              Bez[k].P3[id] = new Point3(Bez_back[k].P3[id].x,
		                       Bez_back[k].P3[id].y,Bez_back[k].P3[id].z);
		      Bez[k].p3[id] = tranfrm(Bez[k].P3[id]);}
		    }
			}
			repaint();
		} 
	else if (evt.target == b5) {
		    morphmode();
			repaint();
		} 
	
		if(evt.target instanceof Checkbox) {
             Checkbox c = (Checkbox)(evt.target);
			 if(c.getLabel().equals("deform")){ 
				 addcirc=0;inc_circ=0;
				  if(findp==1) {
				  }
				  else {findp = 1; n=0;}
				  view_chang=0;} 
			 else if(c.getLabel().equals("add"))
			 { 	 findp = 0; view_chang=1; n=0;
				 addcirc=1;
			 } 
			 else if(c.getLabel().equals("test")){
				 itest = 1-itest;
			 }
		}
		if (evt.target instanceof Choice) {
      //      Choice c = (Choice)(evt.target);
			n = 0; curve1=0;
			switch(cn1.getSelectedIndex()) {				
				case 0 :	Nc = 4;	Ncurve=0;break;
				case 1 :	Nc = 5;	Ncurve=0;break;
				case 2 :	Nc = 6;	Ncurve=0;break;
				case 3 :	Nc = 7;	Ncurve=0;break;
				case 4 :	Nc = 8;	Ncurve=0;break;
				case 5 :	Nc = 9;	Ncurve=0;break;
			}
			switch(cmod.getSelectedIndex()) {
				case 0 :	disp_mod=0;	break;
				case 3 :	fill_crc=1;	break; //Fill circles
				case 4 :	disp_mod=0;	// clear mode
					n = 0;      curve1=0; 
					Ncurve=0;   ncirc=0; findp=0;
					fill_crc=0; disp_anime=0; morph=0; 
					Nmesh=0; mesh_move=0;
					break;
				case 1 :	disp_mod=0;	
				 datafile = "uma1.MET";ncirc=0;
                 scale =0.24; xcenter=40.; ycenter=70.;
		    	 showStatus("Read Data-file Start!");
			     FileInp();
				 inc_circ=1; fill_crc=0; 
			     showStatus("ReadEnd; No. of Metaball data ="+ ncirc);
				 break;
				case 2 :	disp_mod=0;	
				 datafile = "tora2.met";
                 scale =0.24; xcenter=40.; ycenter=50.;
		    	 showStatus("Read Data-file(tora2.met) Start!");
			     FileInp();
				 inc_circ=1; fill_crc=0; 
			     showStatus("ReadEnd; No. of Metaball data ="+ ncirc);
				 repaint();
				 break;
				case 5 :	disp_anime=1; n=0; fill_crc=1;
			      showStatus("animation=" );		
					break; //animation

				case 6 :	disp_mod=2;	break; //Data print
				case 7 :	Nmesh=15;Mmesh=13;mesh_move=0;
			      showStatus("Mesh="+Nmesh );		
					break; //Mesh
            }
			repaint();
		 }
			return true;
	}
/* **************************************************************: */
	public void morphmode() {
		int i,k;

		morph ++;
        System.err.println(" morph="+morph);
        if(morph==1) {
         ncirc_s = ncirc;
		 metaball_s = new Metaball[ncirc+1];
		 circp_s = new Vec2[ncirc+1][5];
		 for( i=1;i<=ncirc_s; i++) {
		 metaball_s[i] = new Metaball(metaball[i].cen.x,
			 metaball[i].cen.y,metaball[i].cen.z,metaball[i].r,metaball[i].v);
		 for( k=0;k<=Ncurve;k++) {
			 circp_s[i][k] = new Vec2(circp[i][k].u, circp[i][k].v);}		 
//		 System.err.println(" i="+i+" u="+circp_s[i][0].u+" v="+circp_s[i][0].v);
		 }
		 for( k=0;k<=Ncurve;k++) {
			Bez_s[k] = new Curve3();
		    Bez_s[k].SetOrg(x_origin,y_origin);
			for(i=0;i<Nc;i++) {
              Bez_s[k].P3[i] = new Point3(Bez[k].P3[i].x,
		                       Bez[k].P3[i].y,Bez[k].P3[i].z);
		      Bez_s[k].p3[i] = tranfrm(Bez_s[k].P3[i]);
			}
			Bez_s[k].copy(Nc);			
		 }		 
         System.err.println(" copy end; ncirc_s="+ncirc_s);
		}
		else if(morph==2) {
         ncirc_d = ncirc;
		 metaball_d = new Metaball[ncirc+1];
		 circp_d = new Vec2[ncirc+1][5];
		 for(i=1;i<=ncirc_d; i++) {
		 metaball_d[i] = new Metaball(metaball[i].cen.x,
			 metaball[i].cen.y,metaball[i].cen.z,metaball[i].r,metaball[i].v);
		 for(k=0;k<=Ncurve;k++) {
			 circp_d[i][k] = new Vec2(circp[i][k].u, circp[i][k].v);}		 
		 }
		 for( k=0;k<=Ncurve;k++) {
			Bez_d[k] = new Curve3();
		    Bez_d[k].SetOrg(x_origin,y_origin);
			for(i=0;i<Nc;i++) {
              Bez_d[k].P3[i] = new Point3(Bez[k].P3[i].x,
		                       Bez[k].P3[i].y,Bez[k].P3[i].z);
		      Bez_d[k].p3[i] = tranfrm(Bez_d[k].P3[i]);
			}
			Bez_d[k].copy(Nc);	
		 }
          System.err.println(" copy end; ncirc_d="+ncirc_d
			+" Ncv="+Ncurve);
		}
}
/* **************************************************************: */
	public void preMorp(Graphics g) {
	if(morph==2) {
	  if(Ncurve==0) taiou();
	  else taiou2(g);
	}
    if(morph > 2) inbtween(g);
}
/* **************************************************************: */
	public Point2  new_point(int i,Vec2 circp_i[][]) {
	 Point3 P0= new Point3(0.,0.,0.);
	 Point3 P1 =new Point3(0.,0.,0.);
	 double a,b,d,x,y;
	 double weight, total_weight;
	 double x1 = 0.;double y1=0.; 

         total_weight =0.; x=0.;y=0.;
		 for (int k=0;k<=Ncurve;k++){
				  double u = circp_i[i][k].u;
				  P0 = (Point3)Bez[k].eval(Nc,u);  
				  P1 = (Point3)Bez[k].eval(Nc,u+0.01);
				  a = -(P1.z-P0.z); b =(P1.y-P0.y);
				  d = Math.sqrt(a*a+b*b); a = a/d; b=b/d; 
				   double v = circp_i[i][k].v;
				   x1 =  P0.y+ a * v;  y1 =  P0.z+ b * v;
				   if(v<0.0) v = -v;
				   weight = 1./(1.+v/10.);  if(v<0.5) weight = 0.;
				   total_weight += weight;
				   x = x+weight*x1; y = y+weight*y1;      
		 }
		return new Point2(x/total_weight,y/total_weight);
}
/* **************************************************************: */
	public void taiou2(Graphics g ) {
		double us,vs,dmin,du,dv,dd;
		int i,j,k,jmin,ir;
		int[] work = new int[ncirc_d+1];

		morph=3;
		g.setColor(Color.white); g.drawString("Waite a moment:",20,65);
        System.err.println(" taiou2="+ncirc_s);
		cores = new int[ncirc_s+ncirc_d+1];
		for( j=1;j<=ncirc_d; j++) work[j]=0;

		 for( i=1;i<=ncirc_s; i++) {
			 Point2 ps = new_point(i,circp_s);
//			 if(i%10==0)
  //      System.err.println(" i="+i+" x="+ps.x+" y="+ps.y+"ncircd="+ncirc_d);
         		g.setColor(Color.green);
				ir=(int)(metaball_s[i].r*half1);
				double x = ps.x+x_origin;
				double y = ps.y+y_origin;
				if(fill_crc == 1) 
					g.fillOval((int)x-ir,(int)y-ir,ir*2,ir*2);
				else g.drawOval((int)x-ir,(int)y-ir,ir+ir,ir+ir);
  			    metaball[i] = new Metaball(x, y,metaball_s[i].cen.z,
					metaball_s[i].r,1.);
                metaball[i].cen.x = x;
                metaball[i].cen.y = y;
           
			 dmin = 1000.; jmin=0;
		   for( j=1;j<ncirc_d; j++) {
			   du = metaball_d[j].cen.x - x;
               dv = metaball_d[j].cen.y - y;
			   dd = Math.sqrt(du*du+dv*dv);
			   if(dd<dmin) { dmin = dd; jmin=j;}
		   }
   //     System.err.println(" i="+i+" j="+jmin+" dmin="+dmin+" x="+metaball_d[jmin].cen.x);
		   cores[i]=jmin;
		   work[jmin]=i;
               g.setColor(Color.yellow);
			   ir=(int)(metaball_d[jmin].r*half1);
				if(fill_crc == 1) 
					g.fillOval((int)metaball_d[jmin].cen.x-ir,
					(int)metaball_d[jmin].cen.y-ir,ir*2,ir*2);
				else g.drawOval((int)metaball_d[jmin].cen.x-ir,
					(int)metaball_d[jmin].cen.y-ir,ir+ir,ir+ir);				
		 }
		   int nadd =0;
		   for( j=1;j<=ncirc_d; j++) {
			   cores[ncirc_s+j]=0;
			   if(work[j]==0){

			  double ud =  metaball_d[j].cen.x;
			  double vd =  metaball_d[j].cen.y;

			  dmin = 1000.; jmin=0;
		      for( i=1;i<=ncirc_s; i++) {
				  /*
			   Point2 ps = new_point(i,circp_s);
			   du = ps.x - ud;   dv = ps.y - vd;
			   */
			   du = metaball[i].cen.x-ud; dv = metaball[i].cen.y-vd;
			   dd = Math.sqrt(du*du+dv*dv);
			   if(dd<dmin) {dmin = dd; jmin=i;}
			  }
		       nadd ++;
		       cores[ncirc_s+j]=-jmin;
			   /*
			   if(nadd%5==0)
             System.err.println(" j="+j+" work="+work[j]+" i="+jmin
			 +" nadd="+nadd+" dmin="+dmin); */
		   }
		   }
		g.setColor(Color.white); g.drawString("Push morph:",20,75);
	}
/* **************************************************************: */
	public void taiou() {
		double us,vs,dmin,du,dv,dd;
		int i,j,k,jmin;
		int[] work = new int[ncirc_d+1];

 		cores = new int[ncirc_s+ncirc_d+1];
		for( j=1;j<=ncirc_d; j++) work[j]=0;
 
		 for( i=1;i<=ncirc_s; i++) {
			 k =0;
			 us = circp_s[i][k].u;
			 vs = circp_s[i][k].v;

			 dmin = 1000.; jmin=0;
		   for( j=1;j<=ncirc_d; j++) {
			   du = circp_d[j][k].u-us;
               dv = circp_d[j][k].v-vs;
			   dd = Math.sqrt(du*du+dv*dv);
			   if(dd<dmin) {
				   dmin = dd; jmin=j;}
		   }
   //     System.err.println(" i="+i+" j="+jmin+" dmin="+dmin);
		   cores[i]=jmin;
		   work[jmin]=i;		   
		 }
		   int nadd =0;
		   for( j=1;j<=ncirc_d; j++) {
			   cores[ncirc_s+j]=0;

			   if(work[j]==0){
			   k =0;
			  double ud = circp_d[j][k].u;
			  double vd = circp_d[j][k].v;

			  dmin = 1000.; jmin=0;
		      for( i=1;i<=ncirc_s; i++) {
			   du = circp_s[i][k].u-ud;
               dv = circp_s[i][k].v-vd;
			   dd = Math.sqrt(du*du+dv*dv);
			   if(dd<dmin) {dmin = dd; jmin=i;}
			  }
		       nadd ++;
		       cores[ncirc_s+j]=-jmin;
       //      System.err.println(" j="+j+" work="+work[j]+" i="+jmin
		//	 +" nadd="+nadd+" dmin="+dmin);
		   }
		   }
	}
/* **************************************************************: */
	public void inbtween(Graphics g) {
		int i,k_step,i1,j1;
		int nstep = 10;
		double t,tt,x,y,z,r,v;

		k_step = morph-3;
		if(k_step>nstep) {morph=3; fill_crc=1-fill_crc; return;}
		
		tt = (double)k_step/(double)nstep;	t = 1.-tt;
		ncirc=0;
		for(i=1;i<=ncirc_s+ncirc_d;i++){
			i1 = i;
			j1 = cores[i];
			if(j1 != 0) { 
			if(j1<0) { j1=i-ncirc_s; i1 = -cores[i];}
	         if(i1>ncirc_s) i1 = ncirc_s;
              if(j1>=ncirc_d) j1 = ncirc_d;
			    x =t*metaball_s[i1].cen.x+tt*metaball_d[j1].cen.x;
				y =t*metaball_s[i1].cen.y+tt*metaball_d[j1].cen.y;
				z =metaball_s[i1].cen.z;
				r =t*metaball_s[i1].r+tt*metaball_d[j1].r;
				ncirc++;
/*				if(ncirc>ncirc_s){
				if(ncirc%20==0)
                System.err.println(" i="+i+" i1="+i1+" j1="+j1+" ncirc="+ncirc);
				}*/
		        metaball[ncirc] = new Metaball(x,y,z,r,1.0);
 //       System.err.println(" i="+i+" j1="+j1+" x="+x+" y="+y+" r="+r);
			}
		}
        System.err.println(" kstep="+k_step+" tt="+tt+" ncirc="+ncirc);
		g.drawString("step="+k_step,20,75);
	}
/* **************************************************************: */
	public void crc_anime() {
      if(disp_anime==1){
			for(int k=0; k<=Ncurve;k++) {
				Bez1[k] = new Curve3();
			for(int id=0;id<Nc;id++) {
              Bez1[k].P3[id] = new Point3(Bez[k].P3[id].x,
		                       Bez[k].P3[id].y,Bez[k].P3[id].z);
		      Bez1[k].p3[id] = tranfrm(Bez1[k].P3[id]);}
		    }
	  }
			  System.err.println(" Disp_anime="+disp_anime);
	  for(int itime=0; itime<120; itime++){
	  double  t= 1.-(double)disp_anime/50.;

	  for(int k=0; k<=Ncurve;k++) {
			for(int id=0;id<Nc;id++) {
              Bez[k].P3[id] = new Point3(t*Bez1[k].P3[id].x+(1.-t)*Bez_back[k].P3[id].x,
		                       t*Bez1[k].P3[id].y+(1.-t)*Bez_back[k].P3[id].y,
							   t*Bez1[k].P3[id].z+(1.-t)*Bez_back[k].P3[id].z);
		      Bez[k].p3[id] = tranfrm(Bez[k].P3[id]);}
		    }
        disp_anime++; idrag=1;
		if(disp_anime>50) {disp_anime=2;idrag=1;}
	  	     showStatus("animation=" +disp_anime);
	   	repaint();}
	}

} // end of "circmor" class