//■二次曲面（quadric surface）
//
//   �@textbox1〜5には次のように設定します。
//
//    　　入力項目名     textBox名  初期値  
//   　　----------------------------------
//    　　拡散反射係数   textBox1　  0.7
//    　　鏡面反射係数   textBox2　  0.4
//    　　反射率　　　   textBox3　  0.4
//    　　環境光係数　   textBox4　  0.2
//    　　屈折率         textBox5　  1.03
//
//   �Alabel5は表面色入力用です。
//   　Clickイベントハンドラにlabel5_Clickを指定します。
//   　また、以下のようにプロパティを設定します。
//   
//
//     AutoSize    False
//     BackColor   Red
//     BorderStyle Fixed3D          
//     Text        (空文字列）
//
//   �Bbotton1〜9にはClickイベントハンドラにそれぞれ
//     <ボタン名>_Clickを指定します。
//
//       ボタン名   Textプロパティ
// 　　------------------------------------
//       button1    楕円面
//       button2    一葉双曲面
//       button3    二葉双曲面
//       button4    楕円錐面
//       button5    楕円放物面
//       button6    双曲放物面
//       button7    放物柱面
//       button8    楕円柱面
//       button9    双曲柱面
//

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using System.Drawing.Drawing2D;


namespace quadricSurface
{
    public partial class Form1 : Form
    {
        private const float 交点なし = 1E8F;
        private Image image;
        private Point org = new Point();
        private TMatrix t, it;
        private Point v0;
        private float dv;
        private float size_x;
        private float size_y;
        private float[,] z_buff;
        private float 拡散反射係数, 鏡面反射係数, 環境光係数;
        private float 反射率, 屈折率;
        Color 物体色;
        public Form1()
        {
            InitializeComponent();
            SetStyle(ControlStyles.DoubleBuffer | ControlStyles.UserPaint |
            ControlStyles.AllPaintingInWmPaint, true);
            image = new Bitmap(500, 400);
            z_buff = new float[image.Width, image.Height];
            Graphics g = Graphics.FromImage(image);
            Brush brush = new SolidBrush(Color.White);
            g.FillRectangle(brush, 0, 0, image.Width, image.Height);
	
		
        }
        protected override void OnPaint(PaintEventArgs e)
        {
            base.OnPaint(e);
            e.Graphics.DrawImage(image, 0, 0);
        }
        public float len(Point a){ return (float)Math.Sqrt(a.X*a.X+a.Y*a.Y+a.Z*a.Z);}
		public Point unit(Point a){ return (len(a)==0)?new Point(1,0,0):a/len(a);}
		public Point unit(float a, float e)
		{
			double PIR=Math.PI/180;
			double aa=a*PIR;
			double ee=e*PIR;
			float X=(float)(Math.Cos(aa)*Math.Cos(ee));
			float Y=(float)(Math.Sin(aa)*Math.Cos(ee));
			float Z=(float)(Math.Sin(ee));
			return new Point(X,Y,Z);
		}
		public Point unit(float a)
		{
			double PIR=Math.PI/180;
			double aa=a*PIR;
			float X=(float)Math.Cos(aa);
			float Y=(float)Math.Sin(aa);
			float Z=0;
			return new Point(X,Y,Z);
		}
		public static TMatrix scale(float scale_x, float scale_y, float scale_z)
		{
			TMatrix m=new TMatrix();
			m.t[0,0]=scale_x; m.t[1,1]=scale_y; m.t[2,2]=scale_z; return m;
		}
		public static TMatrix scale(Point P){ return scale(P.X, P.Y, P.Z);}
		public static TMatrix scale(float S){ return scale(S, S, S);}
		public TMatrix rot_x(float sn, float cs) //X軸を中心に回転
		{
			TMatrix m=new TMatrix();
			m.t[1,1]=cs; m.t[1,2]=-sn; m.t[2,1]=sn; m.t[2,2]=cs; 
			return m;
		}
		public TMatrix rot_y(float sn, float cs) //Y軸を中心に回転
		{
			TMatrix m=new TMatrix();
			m.t[0,0]=cs; m.t[0,2]=sn; m.t[2,0]=-sn; m.t[2,2]=cs; 
			return m;
		}
		public TMatrix rot_z(float sn, float cs) //Z軸を中心に回転
		{
			TMatrix m=new TMatrix();
			m.t[0,0]=cs; m.t[0,1]=-sn; m.t[1,0]=sn; m.t[1,1]=cs; 
			return m;
		}
		public TMatrix rot_x(float r) //X軸を中心に回転
		{
			double PIR=Math.PI/180;
			return rot_x((float)Math.Sin(r*PIR),(float)Math.Cos(r*PIR));
		}
		public TMatrix rot_y(float r) //Y軸を中心に回転
		{
			double PIR=Math.PI/180;
			return rot_y((float)Math.Sin(r*PIR),(float)Math.Cos(r*PIR));
		}
		public TMatrix rot_z(float r) //Z軸を中心に回転
		{
			double PIR=Math.PI/180;
			return rot_z((float)Math.Sin(r*PIR),(float)Math.Cos(r*PIR));
		}
		public TMatrix rotate(Point a, float sn, float cs) //回転変換
		{
			Point P=new Point(a.X, a.Y,0);
			Point axy=unit(P);
			Point axz=rot_z(-axy.Y,axy.X)*unit(a);
			return rot_z(axy.Y, axy.X)*rot_y(-axz.Z,axz.X)*rot_x(sn,cs)
				*rot_y(axz.Z,axz.X)*rot_z(-axy.Y, axy.X);
		}
		public TMatrix rotate(float r) { return rot_z(r);}
		public TMatrix rotate(Point a, float r)
		{
			double PIR=Math.PI/180;
			return rotate(a,(float)Math.Sin(r*PIR),(float)Math.Cos(r*PIR));
		}
		public TMatrix rotate(Point a, Point b)
		{
			float sn=len(unit(a) % unit(b));
			if(sn !=0) return rotate(unit(a%b),sn, unit(a)*unit(b));
			return new TMatrix();
		}
		public TMatrix rotate(Point n) //ベクトル方向に回転
		{
			Point   P  = new Point(1,0,0);
			Point   ss = P%unit(n); if(len(ss)==0) ss= new Point(0,0,1);
			Point   s  = unit(ss);
			TMatrix t1 = rot_x(-s.Y,s.Z);
			TMatrix t3 = rot_x(s.Y,s.Z);
			Point   m  = unit(t1*n);
			TMatrix t2 = rot_z(m.Y, m.X);
			return  t3*t2*t1;
		}
		public TMatrix move(float x, float y, float z)
		{
			TMatrix m=new TMatrix();
			m.t[0,3]=x; m.t[1,3]=y; m.t[2,3]=z;
			return m;
		}
		public TMatrix move(Point P) { return move(P.X, P.Y, P.Z); }
		public TMatrix perspect(float d)
		{
			TMatrix m=new TMatrix();
		    m.t[2,3]=-d; m.t[3,2]=-1/d; return m;
		}
		public TMatrix trans(TMatrix a)
		{
			TMatrix c=new TMatrix();
			for(int i=0; i<4; i++)
				for(int j=0; j<4; j++) c.t[i,j]=a.t[j,i];
			return c;
		}
		public TMatrix inverse(TMatrix a)
		{
			float d=0;
			for(int i=0; i<3; i++)
				d = d + a.t[0,i]*a.t[1, (i+1)%3]*a.t[2, (i+2) % 3]
					- a.t[0,i]*a.t[1, (i+2)%3]*a.t[2, (i+1) % 3];
			TMatrix b=new TMatrix();
			if(d != 0)
				for(int i=0; i<3; i++)
					for(int j=0; j<3; j++)
						b.t[i,j]=(a.t[(j+1)%3,(i+1)%3]*a.t[(j+2)%3, (i+2)%3]
							-a.t[(j+1)%3,(i+2)%3]*a.t[(j+2)%3, (i+1)%3])/d;
			for(int i=0; i<3; i++) b.t[i,3]=-a.t[i,3];
			return b;

		}
		private Surface revolve_z(Line a, int n,float b, float e)//ｚ軸回転によるポリゴンモデル作成
		{
			Surface f=new Surface((a.n-1)*n);
			TMatrix m0,m1;
			for(int j=0,t=0;j<n; j++)
			{
				m0=rot_z(b+j*(e-b)/n);
				m1=rot_z(b+(j+1)*(e-b)/n);
				for(int i=0;i<a.n-1;i++, t++)
					f.p[t]=new Polygon(m0*a.p[i],m0*a.p[i+1],m1*a.p[i+1],m1*a.p[i]);
			}
			return f;
		}
		private Surface revolve_z(Line a, int n,float b)//ｚ軸回転によるポリゴンモデル作成
		{
			return revolve_z(a, n, b, 360);
		}
		private Surface revolve_z(Line a, int n)//ｚ軸回転によるポリゴンモデル作成
		{
			return revolve_z(a, n, 0);
		}
		public TMatrix rot_xz(Point a)
		{
			Point ua=unit(a);
			Point aa=scale(1,1,0)*ua;
		    float lenAA=len(aa);
			if(lenAA==0) return rot_x(90);
			TMatrix mx=rot_x((new Point(0,0,1))*ua, lenAA);
			TMatrix mz=rot_z(-aa.X/lenAA, aa.Y/lenAA);
			return mz*mx;
		}
		public Surface sweep_xz(Line a, Line b)
		{
			//Point[] p=new Point[(a.n)*(b.n)];
			Line[]  f=new Line[b.n];
			//for(int j=0; j<b.n; j++) f[j].p=p[j*a.n];
			if(len(b.p[0]-b.p[b.n-1])>0.1)
			{
				f[0]    =move(b.p[0])    *rot_xz(b.p[1]-b.p[0])        *a;
				f[b.n-1]=move(b.p[b.n-1])*rot_xz(b.p[b.n-1]-b.p[b.n-2])*a;
			}
			else
			{
				f[0]    =move(b.p[0])*rot_xz(b.p[1]-b.p[b.n-4])*a;
				f[b.n-1]=move(b.p[0])*rot_xz(b.p[3]-b.p[b.n-2])*a;
			}
			for(int i=1;i<b.n-1; i++)
				f[i]=move(b.p[i])*rot_xz(b.p[i+1]-b.p[i-1])*a;
			Surface s=new Surface((a.n-1)*(b.n-1));
			for(int i=0,t=0;i<b.n-1; i++)
				for(int j=0;j<a.n-1;j++,t++)
					s.p[t]=new Polygon(f[i].p[j],f[i].p[j+1],f[i+1].p[j+1],f[i+1].p[j]);
            return s;

		}
		private void init2D(int x, int y)
		{
			//SetStyle(ControlStyles.DoubleBuffer | ControlStyles.UserPaint |
			//	ControlStyles.AllPaintingInWmPaint, true);
			try
			{
				image=new Bitmap(x, y);
				Graphics g =Graphics.FromImage(image);
				Brush brush = new SolidBrush(Color.White);
				g.FillRectangle(brush, 0, 0, image.Width,image.Height);
				size_x=x*2; 
				size_y=y*2;
				org.X=(float)(size_x/2-0.5);
				org.Y=(float)(size_y/2-1.5);
				org.Z=0;
				setview(30,30,1000);
				this.Invalidate();
			}
			catch ( Exception err)
			{
				MessageBox.Show(err.Message);
			}
		}
		private void origin(int x, int y)
		{
			org=new Point((float)x+0.5F , (float)y-0.5F,0);
		}
		private void debugTmat(string Name,TMatrix t)
		{
			string S="*** "+  Name + " ***\r\n"; 
			for(int i=0;i<4;i++)
			{   S=S+"\r\n";
				for(int j=0;j<4;j++)
				{
					S=S +"\t"+t.t[i,j];
				}
			}
			MessageBox.Show(S);
		}
		private void debugPoint(string Name,Point P)
		{
			string S="*** "+  Name + " ***\r\n"; 
			  S=S+" X = " + P.X +"\t Y = "+ P.Y+"\t Z = "+ P.Z ;
				
			MessageBox.Show(S);
		}

		private void setview(float a,float e, float d)//始点方向設定
		{
			Point P1=new Point(0,0,1);
			dv=d;
			t=perspect(d)*rot_x(e-90)*rot_z(-90-a);//	debugTmat("透視変換",t);
			it=rot_z(90+a)*rot_x(90-e);
			v0=it*P1*d;
		}
		private void setview(Point p)//始点方向設定
		{
			Point P1=new Point(0,0,1);
			dv=len(p);
			Point wa=unit(p); v0=p;
			if((wa.X==0) &&(wa.Y==0)){ t=perspect(dv);return;}
            float wd=(float)Math.Sqrt(wa.X*wa.X+wa.Y*wa.Y);
            t=perspect(dv)*rot_x(-wd,wa.Z)*rot_z(-wa.X/wd,-wa.Y/wd);//	debugTmat("透視変換",t);
			it=rot_z(wa.X/wd,-wa.Y/wd)*rot_x(wd,wa.Z);
		}
		private Point screen(Point p)
		{
			Point a=t*p;
			//debugTmat("screen t",t);
			//debugPoint("screen p",p);
			//debugPoint("screen a",a);
			
			float s=t.t[3,0]*p.X+t.t[3,1]*p.Y+t.t[3,2]*p.Z+t.t[3,3];
			//MessageBox.Show(" s=" + s.ToString());
			return new Point(org.X+a.X/s,org.Y-a.Y/s,a.Z);
		}
		
		private Point world(Point p) //View座標をWorld座標に変換
		{
			float sz=-p.Z/dv;
			return it*(new Point((p.X-org.X)*sz,(-p.Y+org.Y)*sz,p.Z+dv));
		}
		
		

		private void drawAxis()
		{
			Graphics g = Graphics.FromImage(image);　 // 表示先の指定 
			Pen pen= new Pen(Color.Black);　　　　　　// 座標軸の描画
			pen.EndCap=LineCap.ArrowAnchor;
			Font font =new Font(Font.Name,10);
			Brush brush=new SolidBrush(Color.Black);

  　　      drawAxis(g, pen,"X",font,brush,new Point((float)(-size_x/2),0,0), new Point((float)size_x/2,0,0));
  　　      drawAxis(g, pen,"Y",font,brush,new Point(0,(float)(-size_x/2),0), new Point(0,(float)(size_x/2*0.9),0));
  　　      drawAxis(g, pen,"Z",font,brush,new Point(0,0,(float)(-size_y/2)), new Point(0,0,(float)(size_y/2)));
		}
		private void drawAxis(Graphics g, Pen pen, string AX, Font font,Brush brush,Point P1, Point P2)
		{

  	　　　  Point PP1 = screen(P1);
			Point PP2 = screen(P2);
	
			//標準座標からスクリーン座標への変換
　　　　　　g.DrawLine(pen, PP1.X, PP1.Y, PP2.X, PP2.Y);
			g.DrawString(AX,font,brush,PP2.X,PP2.Y);

		}
		private void drawLine(Graphics g, Pen pen, Point P1, Point P2)
		{
	　　　  Point PP1 = screen(P1);
			Point PP2 = screen(P2);
	
			//標準座標からスクリーン座標への変換
　　　　　　g.DrawLine(pen, PP1.X, PP1.Y, PP2.X, PP2.Y);
		}
		private void drawLine(Point P1, Point P2,Color C)
		{
			Graphics g = Graphics.FromImage(image);　 // 表示先の指定 
			Pen pen= new Pen(C);　　　　　　// 座標軸の描画
			pen.EndCap=LineCap.NoAnchor;
			Font font =new Font(Font.Name,20);
			Brush brush=new SolidBrush(Color.Black);
            
	　　　  Point PP1 = screen(P1);
			Point PP2 = screen(P2);
	
			//標準座標からスクリーン座標への変換
　　　　　　g.DrawLine(pen, PP1.X, PP1.Y, PP2.X, PP2.Y);
		}
		private void drawLineDot(Point P1, Point P2,Color C)
		 {
			 Graphics g = Graphics.FromImage(image);　 // 表示先の指定 
			 Pen pen= new Pen(C);　　　　　　// 座標軸の描画
			 pen.DashStyle=DashStyle.Dot;
			 pen.EndCap=LineCap.NoAnchor;
			 Font font =new Font(Font.Name,20);
			 Brush brush=new SolidBrush(Color.Black);
            
	　　　  Point PP1 = screen(P1);
			 Point PP2 = screen(P2);
	
			 //標準座標からスクリーン座標への変換
　　　　　　g.DrawLine(pen, PP1.X, PP1.Y, PP2.X, PP2.Y);
		 }
		private void drawLine(Line LDT,Color C)
		{
			Point P1=LDT.p[0];
			for(int i=1;i<LDT.n;i++)
			{
				Point P2=LDT.p[i];
				drawLine(P1, P2, C);
				P1=P2;
			}
		}
		private void drawLine(Polygon LDT,Color C)
		{
			Point P1=LDT.p[0];
			for(int i=1;i<LDT.n;i++)
			{
				Point P2=LDT.p[i];
				drawLine(P1, P2, C);
				P1=P2;
			}
		}
		private void drawLine(Surface f,Color C)
		{
			for(int i=0; i<f.n; i++) drawLine(f.p[i],C);
		}
		
		private void paintPolygon(Polygon p,Color C)
		{
			System.Drawing.PointF []q= new System.Drawing.PointF[p.n];
			Point a;
			for(int i=0;i<p.n;i++){ a=screen(p.p[i]);q[i].X=a.X;q[i].Y=a.Y;}
			Graphics g = Graphics.FromImage(image);　 // 表示先の指定 
			
			Brush brush = new SolidBrush(C);
			g.FillPolygon(brush,q);
		}
		private void drawPoint(float x,float y,Color C)
		{
			Graphics g = Graphics.FromImage(image);　 // 表示先の指定 
			Brush brush = new SolidBrush(C);
			g.FillRectangle(brush, x, y, 1F, 1F);
		
		}
		private void drawPoint(Point p,Color C)
		{
			drawPoint(p.X,p.Y,C);
		
		}
		public Line circle_xz(int n, float b, float e)
		{
			Line c=new Line(n+1);
			double PIR=Math.PI/180;
			double TH=b*PIR;double DTH=PIR*(e-b)/n;
			for(int i=0; i<n+1;i++,TH+=DTH) c.p[i]=new Point((float)Math.Sin(TH),0,(float)Math.Cos(TH));
			return c;
		}
		public Line circle_xz(int n){return circle_xz(n,0,360);}
		public Line circle_yz(int n, float b, float e)
		{
			Line c=new Line(n+1);
			double PIR=Math.PI/180;
			double TH=b*PIR;double DTH=PIR*(e-b)/n;
			for(int i=0; i<n+1;i++,TH+=DTH) c.p[i]=new Point(0,(float)Math.Cos(TH),(float)Math.Sin(TH));
			return c;
		}
		public Line circle_yz(int n){return circle_yz(n,0,360);}
		public Line circle_xy(int n, float b, float e)
		{
			Line c=new Line(n+1);
			double PIR=Math.PI/180;
			double TH=b*PIR;double DTH=PIR*(e-b)/n;
			for(int i=0; i<n+1;i++,TH+=DTH) c.p[i]=new Point((float)Math.Cos(TH),(float)Math.Sin(TH),0);
			return c;
		}
		public Line circle_xy(int n){return circle_xy(n,0,360);}

		//private Point[] 線分変換(TMatrix TM, Point[] LD)
		//{
		//	int N=LD.Length;
		//	Point [] newLine= new Point[N];
		//	for(int i=0;i<N;i++)
		//	{
		//		newLine[i]=TM*LD[i];                 // 点変換(TM,LD[i]);
		//	}
		//	return newLine;
		//}
		//反射光の計算 V:視線ベクトル  L:照明光ベクトル N：表面の法線ベクトル
		private Color shading(Point V, Point L, Point N, Color CU)
		{
            double kd=拡散反射係数, ks=鏡面反射係数, ke=環境光係数;
			Point RV=L-2*(L*N)*N;            //反射方向ベクトル
			
			double D1=(-N*L)/(len(N)*len(L));　 // 拡散および環境光
			if(D1<0)D1=0;
			D1=kd*D1+ke;　                         
			
			double D2=ks*((-RV*V))/(len(RV)*len(V));　　　   // 反射
			if(D2<0)D2=0;
			int R=(int)(D1*(double)(CU.R)+D2*255);if(R>255)R=255;
			int G=(int)(D1*(double)(CU.G)+D2*255);if(G>255)G=255;
			int B=(int)(D1*(double)(CU.B)+D2*255);if(B>255)B=255;
            return Color.FromArgb(R,G,B);
		}
		private Color mix(Color C1,double d1,Color C2,double d2)
		{
			int R=(int)(d1*(double)(C1.R)+d2*(double)(C2.R));if(R>255)R=255;
			int G=(int)(d1*(double)(C1.G)+d2*(double)(C2.G));if(G>255)G=255;
			int B=(int)(d1*(double)(C1.B)+d2*(double)(C2.B));if(B>255)B=255;
			return Color.FromArgb(R,G,B);
		}
		private Color mix(Color C1,double d1,Color C2,double d2,Color C3,double d3)
		{
			int R=(int)(d1*(double)(C1.R)+d2*(double)(C2.R)+d3*(double)(C3.R));if(R>255)R=255;
			int G=(int)(d1*(double)(C1.G)+d2*(double)(C2.G)+d3*(double)(C3.G));if(G>255)G=255;
			int B=(int)(d1*(double)(C1.B)+d2*(double)(C2.B)+d3*(double)(C3.B));if(B>255)B=255;
			return Color.FromArgb(R,G,B);
		}
		private Color shading(Point V, Point L, Point N, Color CU, double s)
		{
			double kd=拡散反射係数, ks=鏡面反射係数, ke=環境光係数;
			Point RV=L-2*(L*N)*N;            //反射方向ベクトル
			
			double D1=(-N*L)/(len(N)*len(L));　 // 拡散および環境光
			if(D1<0)D1=0;
			D1=kd*s*D1+ke;　                         
			
			double D2=s*ks*((-RV*V))/(len(RV)*len(V));　　　   // 反射
			if(D2<0)D2=0;
			int R=(int)(D1*(double)(CU.R)+D2*255);if(R>255)R=255;
			int G=(int)(D1*(double)(CU.G)+D2*255);if(G>255)G=255;
			int B=(int)(D1*(double)(CU.B)+D2*255);if(B>255)B=255;
			return Color.FromArgb(R,G,B);
		}
		private Point 補間(Point p1,Point p2,int a1,int a2)//p1,p2を補間
		{
			return ((a1+a2<1) ?p1: ((a2)*p1+a1*p2)/(a1+a2));
		}
		private Color 色補間(Color c1,Color c2,int a1,int a2)//p1,p2を補間
		{
			int aa=a1+a2;
			float b1=(float)a1, b2=(float)a2;
			if(aa<1) return c1;
			float R1=(float)c1.R, G1=(float)c1.G, B1=(float)c1.B;
			float R2=(float)c2.R, G2=(float)c2.G, B2=(float)c2.B;
            
			return Color.FromArgb((int)((b2*R1+b1*R2)/aa),(int)((b2*G1+b1*G2)/aa),(int)((b2*B1+b1*B2)/aa));
		}
			
		public Csgm ball(Point s)
		{
			Point s2=new Point((s.X==0)?0:1/(s.X*s.X),(s.Y==0)?0:1/(s.Y*s.Y),(s.Z==0)?0:1/(s.Z*s.Z));
			//MessageBox.Show("ball X=" +s2.X +" Y=" +s2.Y + " Z=" +s2.Z);
			return new Csgm(s2, new Point(0,0,0),-1);
		}
		public Csgm ball    (float r)          { return ball(new Point(r,r,r));}
		private void 直線描画(Hyper Hy,Point p0, Point p1,Color c)
		{
			Point a0=screen(p0), a1=screen(p1), P=new Point(0,0,0);
			int xs=(a1.X-a0.X>0)? 1: -1;
			int ys=(a1.Y-a0.Y>0)? 1: -1;
			Bitmap w=new Bitmap(image);
			double sx=(double)w.Width/(double)size_x;
			double sy=(double)w.Height/(double)size_y;


			if(Math.Abs((double)(a0.X-a1.X))>Math.Abs((double)(a0.Y-a1.Y)))
			{
				
				int xst=(int)(a0.X+0.5F), xed=(int)(a1.X+0.5F);
				float DD=a1.X-a0.X;
				Point B=new Point(0,0,0);
				for(int x=xst; x!=xed;x+=xs)
				{
					B=a0+(a1-a0)*(x-a0.X)/DD;
					Point A=world(B);
					Point V=unit(A-v0), N=new Point(0,0,0);
					Ray VR=new Ray(v0,V); //視線
					if(Hy.hit(VR,ref P,ref N)>len(A-v0)) drawPoint(x, B.Y, c);
					else drawPoint(x, B.Y, mix(c,0.5, w.GetPixel((int)(sx*(double)x),(int)(sy*(double)B.Y)), 1));
				}
			}
			else
			{
				int yst=(int)(a0.Y+0.5F), yed=(int)(a1.Y+0.5F);
				float DD=a1.Y-a0.Y;
				Point B=new Point(0,0,0);
				
				for(int y=yst; y!=yed;y+=ys)
				{
					B=a0+(a1-a0)*(y-a0.Y)/DD;
					Point A=world(B);
					Point V=unit(A-v0), N=new Point(0,0,0);
					Ray VR=new Ray(v0,V); //視線
					if(Hy.hit(VR, ref P, ref N)>len(A-v0)) drawPoint((int)B.X, y, c);
					else drawPoint(B.X, y, mix(c, 0.5, w.GetPixel((int)(sx*(double)B.X),(int)(sy*(double)y)), 1));
				}
			}

		}
		private void 画面クリア(Color c)
		{
			SetStyle(ControlStyles.DoubleBuffer | ControlStyles.UserPaint |
				ControlStyles.AllPaintingInWmPaint, true);
			Graphics g =Graphics.FromImage(image);
			Brush brush = new SolidBrush(c);
			g.FillRectangle(brush, 0, 0, image.Width,image.Height);
		}

		private void 曲面表示(Point S2,Point S1,float S0)
		{
			Hyper Hy=new Hyper(S2,S1,S0);
			Color C=物体色,	bc=Color.FromArgb(0,0,0);
			
			init2D(500,400); origin(250,200);
			画面クリア(bc);
			setview(60,40,3000);
			Point V0 = v0;                     //視点
			Point L0 = new Point(500,500,500);   //照明光
			Point V  = unit(new Point(0,0,1)-V0), N=new Point(0,0,0), P=new Point(0,0,0);
			Ray   VR = new Ray(V0,V);
			if(Hy.hit(VR, ref P, ref N)>=交点なし)
			{
				if((2*Hy.S2.Z + Hy.S1.Z)>0)
				{ Hy.S2=-Hy.S2; Hy.S1=-Hy.S1; Hy.S0=-Hy.S0;}
				else if(((2*Hy.S2.Z + Hy.S1.Z)==0)&&((2*Hy.S2.Y + Hy.S1.Y)>0))
				{ Hy.S2=-Hy.S2; Hy.S1=-Hy.S1; Hy.S0=-Hy.S0;}
				else if(((2*Hy.S2.Z + Hy.S1.Z)==0)&&((2*Hy.S2.Y + Hy.S1.Y)==0)&&((2*Hy.S2.X + Hy.S1.X)>0))
				{ Hy.S2=-Hy.S2; Hy.S1=-Hy.S1; Hy.S0=-Hy.S0;}
			}
			if(Hy.hit(VR, ref P, ref N)<交点なし)//視線が２次曲面と交わるかどうかの判定
			{
				if(N*V>0) { Hy.S2=-Hy.S2; Hy.S1=-Hy.S1; Hy.S0=-Hy.S0;}
			}
			Bitmap w=new Bitmap(image);
			double sx=(double)w.Width/(double)size_x;
			double sy=(double)w.Height/(double)size_y;
			
			for(int x=0; x<size_x; x++)
				for(int y=0; y<size_y; y++)
				{
					V  = unit(world(new Point(x,y,-dv))-V0);  //視線ベクトル
					VR = new Ray(V0,V);

					if(Hy.hit(VR,ref P, ref N)<交点なし)   //視線が２次曲面と交わるとき
					{
						Point  L = unit(P-L0);             //照明光方向ベクトル

						Color CC=mix(shading(V,L,N,C),1,w.GetPixel((int)(sx*(double)x),(int)(sy*(double)y)),1);
						drawPoint(x,y, CC);
						if(Hy.hit(new Ray(P+V,V), ref P, ref N)<交点なし)
						{
							L = unit(P-L0);             //照明光方向ベクトル
							Color cw=w.GetPixel((int)(sx*(double)x),(int)(sy*(double)y));
							Color cs=shading(V,L,N,C,0.7);
							Color cc=mix(cs,1,cw,1);
							drawPoint(x,y, mix(cs,1,cw,1));
						}
				
					}
				}
            Graphics g = Graphics.FromImage(image);　 // 表示先の指定 
			Font font =new Font(Font.Name,10);
			Brush brush=new SolidBrush(Color.White);

			Point PPX=new Point(180,  0,  0);
			Point PPY=new Point(  0,180,  0);
			Point PPZ=new Point(  0,  0,180);
            Point PP;
  　　      //drawAxis();

			直線描画(Hy,-PPX,PPX,Color.Cyan); PP = screen(PPX);	g.DrawString("X",font,brush,PP.X,PP.Y);
			直線描画(Hy,-PPY,PPY,Color.Cyan); PP = screen(PPY);	g.DrawString("Y",font,brush,PP.X,PP.Y);
			直線描画(Hy,-PPZ,PPZ,Color.Cyan); PP = screen(PPZ);	g.DrawString("Z",font,brush,PP.X,PP.Y);
		
		}

		private void 楕円面()
		{
			float a=40,b=40,c=50,s0;
			Point S2=new Point(1/(a*a), 1/(b*b), 1/(c*c)),S1=new Point(0,0,0); s0=-1;
			曲面表示(S2,S1,s0);
		}
		private void 一葉双曲面()
		{
			float a=40,b=40,c=50,s0;
			Point S2=new Point(1/(a*a), 1/(b*b), -1/(c*c)),S1=new Point(0,0,0); s0=-1;
			曲面表示(S2,S1,s0);
		}

		private void 二葉双曲面()
		{
			float a=40,b=40,c=50,s0;
			Point S2=new Point(-1/(a*a), -1/(b*b), 1/(c*c)),S1=new Point(0,0,0); s0=-1;
			曲面表示(S2,S1,s0);
		}
		private void 楕円錘面()
		{
			float a=40,b=40,c=50,s0;
			Point S2=new Point(1/(a*a), 1/(b*b), -1/(c*c)),S1=new Point(0,0,0); s0=0;
			曲面表示(S2,S1,s0);
		}
		private void 双曲放物面()
		{
			float a=40,b=40,c=50,s0;
			Point S2=new Point(1/(a*a), -1/(b*b), 0),S1=new Point(0,0,-1/c); s0=0;
			曲面表示(S2,S1,s0);
		}
		private void 楕円放物面()
		{
			float a=40,b=40,c=50,s0;
			Point S2=new Point(1/(a*a), 1/(b*b), 0),S1=new Point(0,0,-1/c); s0=0;
			曲面表示(S2,S1,s0);
		}
		
		private void 放物柱面()
		{
			float a=40,b=40,c=50,s0;
			Point S2=new Point(1/(a*a), 0, 0),S1=new Point(0,1/b,0); s0=0;
			曲面表示(S2,S1,s0);
		}
		private void 楕円柱面()
		{
			float a=40,b=40,c=50,s0;
			Point S2=new Point(1/(a*a), 1/(b*b), 0),S1=new Point(0,0,0); s0=-1;
			曲面表示(S2,S1,s0);
		}
		private void	双曲柱面()
		{
			float a=40,b=40,c=50,s0;
			Point S2=new Point(1/(a*a), -1/(b*b), 0),S1=new Point(0,0,0); s0=-1;
			曲面表示(S2,S1,s0);
		}
		private void button1_Click(object sender, System.EventArgs e)
		{
			拡散反射係数=float.Parse(textBox1.Text);
			鏡面反射係数=float.Parse(textBox2.Text);
			環境光係数  =float.Parse(textBox3.Text);
			反射率      =float.Parse(textBox4.Text);
			屈折率      =float.Parse(textBox5.Text);
			物体色      =label5.BackColor;
			楕円面();
			this.Invalidate();

		}
		private void button2_Click(object sender, System.EventArgs e)
		{
			拡散反射係数=float.Parse(textBox1.Text);
			鏡面反射係数=float.Parse(textBox2.Text);
			環境光係数  =float.Parse(textBox3.Text);
			反射率      =float.Parse(textBox4.Text);
			屈折率      =float.Parse(textBox5.Text);
			物体色      =label5.BackColor;

		　　一葉双曲面();
			this.Invalidate();

		}

		private void button3_Click(object sender, System.EventArgs e)
		{
			拡散反射係数=float.Parse(textBox1.Text);
			鏡面反射係数=float.Parse(textBox2.Text);
			環境光係数  =float.Parse(textBox3.Text);
			反射率      =float.Parse(textBox4.Text);
			屈折率      =float.Parse(textBox5.Text);
			物体色      =label5.BackColor;
			
            二葉双曲面();
			this.Invalidate();

		}
		private void button4_Click(object sender, System.EventArgs e)
		{
			拡散反射係数=float.Parse(textBox1.Text);
			鏡面反射係数=float.Parse(textBox2.Text);
			環境光係数  =float.Parse(textBox3.Text);
			反射率      =float.Parse(textBox4.Text);
			屈折率      =float.Parse(textBox5.Text);
			物体色      =label5.BackColor;
			
			楕円錘面();
			this.Invalidate();

		
		}
		private void button5_Click(object sender, System.EventArgs e)
		{
			拡散反射係数=float.Parse(textBox1.Text);
			鏡面反射係数=float.Parse(textBox2.Text);
			環境光係数  =float.Parse(textBox3.Text);
			反射率      =float.Parse(textBox4.Text);
			屈折率      =float.Parse(textBox5.Text);
			物体色      =label5.BackColor;
			
			楕円放物面();
			this.Invalidate();
		
		}
		private void button6_Click(object sender, System.EventArgs e)
		{
			拡散反射係数=float.Parse(textBox1.Text);
			鏡面反射係数=float.Parse(textBox2.Text);
			環境光係数  =float.Parse(textBox3.Text);
			反射率      =float.Parse(textBox4.Text);
			屈折率      =float.Parse(textBox5.Text);
			物体色      =label5.BackColor;
			
			双曲放物面();
			this.Invalidate();
		}
		private void label5_Click(object sender, System.EventArgs e)
		{
			ColorDialog MyDialog = new ColorDialog();
			MyDialog.AllowFullOpen = false ;
			MyDialog.ShowHelp = true ;
			MyDialog.Color = label5.BackColor ;
    		if (MyDialog.ShowDialog() == DialogResult.OK)
				label5.BackColor =  MyDialog.Color;

		}

		private void button7_Click(object sender, System.EventArgs e)
		{
			拡散反射係数=float.Parse(textBox1.Text);
			鏡面反射係数=float.Parse(textBox2.Text);
			環境光係数  =float.Parse(textBox3.Text);
			反射率      =float.Parse(textBox4.Text);
			屈折率      =float.Parse(textBox5.Text);
			物体色      =label5.BackColor;
			
			放物柱面();
			this.Invalidate();
		}

		private void button8_Click(object sender, System.EventArgs e)
		{
			拡散反射係数=float.Parse(textBox1.Text);
			鏡面反射係数=float.Parse(textBox2.Text);
			環境光係数  =float.Parse(textBox3.Text);
			反射率      =float.Parse(textBox4.Text);
			屈折率      =float.Parse(textBox5.Text);
			物体色      =label5.BackColor;
			
			楕円柱面();
			this.Invalidate();
		}

		private void button9_Click(object sender, System.EventArgs e)
		{
			拡散反射係数=float.Parse(textBox1.Text);
			鏡面反射係数=float.Parse(textBox2.Text);
			環境光係数  =float.Parse(textBox3.Text);
			反射率      =float.Parse(textBox4.Text);
			屈折率      =float.Parse(textBox5.Text);
			物体色      =label5.BackColor;
			
			双曲柱面();
			this.Invalidate();
		}
	}

	public class Point
	{
		public float X,Y,Z;
		public Point(float IX,float IY, float IZ){X=IX; Y=IY;Z=IZ;}
		public Point(float IX,float IY)          {X=IX; Y=IY;Z=1; }
		public Point() {X=0;Y=0;Z=0;}
		public static Point operator +(Point a,Point b){return new Point(a.X+b.X, a.Y+b.Y, a.Z+b.Z);}
		public static Point operator -(Point a,Point b){return new Point(a.X-b.X, a.Y-b.Y, a.Z-b.Z);}
		public static Point operator -(Point a)        {return new Point(   -a.X,    -a.Y,    -a.Z);}
		public static Point operator *(float a,Point b){return new Point(  a*b.X,   a*b.Y,   a*b.Z);}
		public static Point operator *(Point a,float b){return new Point(  a.X*b,   a.Y*b,   a.Z*b);}
		public static float operator *(Point a,Point b){return a.X*b.X+ a.Y*b.Y+ a.Z*b.Z;}    //内積

		public static Point operator /(Point a,float b){return new Point(  a.X/b,   a.Y/b,   a.Z/b);}
		
		public static Point operator %(Point a,Point b)                                          //外積 
		{	return new Point(a.Y*b.Z-a.Z*b.Y, a.Z*b.X-a.X*b.Z, a.X*b.Y-a.Y*b.X);}
		public static bool  operator ==(Point a, Point b){ return a.X==b.X && a.Y==b.Y && a.Z==b.Z;}
		public static bool  operator !=(Point a, Point b){ return !(a==b) ;}
		public float len(Point a){ return (float)Math.Sqrt(a.X*a.X+a.Y*a.Y+a.Z*a.Z);}
		public Point unit(Point a){ return (len(a)==0)?new Point(1,0,0):a/len(a);}
		public Point unit(float a, float e)
		{
			double PIR=Math.PI/180;
			double aa=a*PIR;
			double ee=e*PIR;
			float X=(float)(Math.Cos(aa)*Math.Cos(ee));
			float Y=(float)(Math.Sin(aa)*Math.Cos(ee));
			float Z=(float)(Math.Sin(ee));
			return new Point(X,Y,Z);
		}
		public Point unit(float a)
		{
			double PIR=Math.PI/180;
			double aa=a*PIR;
			float X=(float)Math.Cos(aa);
			float Y=(float)Math.Sin(aa);
			float Z=0;
			return new Point(X,Y,Z);
		}
	}
	public class Line : Point 
	{
		
		public int n;
		public Point[] p;
		public Line(int 頂点数){n=頂点数;if(n>0) p=new Point[n];}
		public Line(Point P0, Point P1)          { n=2; p=new Point[n]; p[0]=P0; p[1]=P1;}
		public Line(Point P0, Point P1, Point P2){ n=3; p=new Point[n]; p[0]=P0; p[1]=P1; p[2]=P2;}
		public Line(Point P0,Point P1,Point P2,Point P3)
		{ n=4; p=new Point[n];p[0]=P0;p[1]=P1;p[2]=P2;p[3]=P3;}
		public Line(Point P0,Point P1,Point P2,Point P3, Point P4)
		{ n=5; p=new Point[n];p[0]=P0;p[1]=P1;p[2]=P2;p[3]=P3;p[4]=P4;}
		public Line(Point P0,Point P1,Point P2,Point P3,Point P4,Point P5)
		{ n=6; p=new Point[n];p[0]=P0;p[1]=P1;p[2]=P2;p[3]=P3;p[4]=P4;p[5]=P5;}
		public Line(Point P0,Point P1,Point P2,Point P3,Point P4,Point P5,Point P6)
		{ n=7; p=new Point[n];p[0]=P0;p[1]=P1;p[2]=P2;p[3]=P3;p[4]=P4;p[5]=P5;p[6]=P6;}
		public Line(Point P0,Point P1,Point P2,Point P3,Point P4,Point P5,Point P6,Point P7)
		{ n=8; p=new Point[n];p[0]=P0;p[1]=P1;p[2]=P2;p[3]=P3;p[4]=P4;p[5]=P5;p[6]=P6;p[7]=P7;}
		public Line(Point P0,Point P1,Point P2,Point P3,Point P4,Point P5,Point P6,Point P7,Point P8)
		{ n=9; p=new Point[n];p[0]=P0;p[1]=P1;p[2]=P2;p[3]=P3;p[4]=P4;p[5]=P5;p[6]=P6;p[7]=P7;p[8]=P8;}
		public Line(Point P0,Point P1,Point P2,Point P3,Point P4,Point P5,Point P6,Point P7,Point P8,Point P9)
		{ n=10;p=new Point[n];p[0]=P0;p[1]=P1;p[2]=P2;p[3]=P3;p[4]=P4;p[5]=P5;p[6]=P6;p[7]=P7;p[8]=P8;p[9]=P9;}
		public Line(Line b){ n=b.n;if(n>0) p=new Point[n];for(int i=0;i<n;i++)p[i]=b.p[i];}
        		
		public static Line operator *(float a,Line b)
		{
			Line C=new Line(b);
			for(int i=0;i<C.n;i++) {C.p[i].X *=a;C.p[i].Y *=a;C.p[i].Z *=a;}
			return C;
		}
		public float len(Line a)
		{
			float LN=0; 
			for(int i=0;i<a.n-1;i++) 
			{
				LN += len(a.p[i]-a.p[i+1]);
			}
			return LN;
		}

	}
	public class Polygon : Point 
	{
		public int n;
		public Point[] p;
		public Polygon(int 頂点数){n=頂点数+1;if(n>0) p=new Point[n];}

		public Polygon(Point P0, Point P1){ n=3; p=new Point[n]; p[0]=P0; p[1]=P1; p[2]=P0;}
		public Polygon(Point P0, Point P1, Point P2){ n=4; p=new Point[n]; p[0]=P0; p[1]=P1; p[2]=P2;p[3]=P0;}
		public Polygon(Point P0,Point P1,Point P2,Point P3)
		{ n=5; p=new Point[n];p[0]=P0;p[1]=P1;p[2]=P2;p[3]=P3;p[4]=P0;}
		public Polygon(Point P0,Point P1,Point P2,Point P3, Point P4)
		{ n=6; p=new Point[n];p[0]=P0;p[1]=P1;p[2]=P2;p[3]=P3;p[4]=P4;p[5]=P0;}
		public Polygon(Point P0,Point P1,Point P2,Point P3,Point P4,Point P5)
		{ n=7; p=new Point[n];p[0]=P0;p[1]=P1;p[2]=P2;p[3]=P3;p[4]=P4;p[5]=P5;p[6]=P0;}
		public Polygon(Point P0,Point P1,Point P2,Point P3,Point P4,Point P5,Point P6)
		{ n=8; p=new Point[n];p[0]=P0;p[1]=P1;p[2]=P2;p[3]=P3;p[4]=P4;p[5]=P5;p[6]=P6;p[7]=P0;}
		public Polygon(Point P0,Point P1,Point P2,Point P3,Point P4,Point P5,Point P6,Point P7)
		{ n=9; p=new Point[n];p[0]=P0;p[1]=P1;p[2]=P2;p[3]=P3;p[4]=P4;p[5]=P5;p[6]=P6;p[7]=P7;p[8]=P0;}
		public Polygon(Point P0,Point P1,Point P2,Point P3,Point P4,Point P5,Point P6,Point P7,Point P8)
		{ n=10; p=new Point[n];p[0]=P0;p[1]=P1;p[2]=P2;p[3]=P3;p[4]=P4;p[5]=P5;p[6]=P6;p[7]=P7;p[8]=P8;p[9]=P0;}
		public Polygon(Point P0,Point P1,Point P2,Point P3,Point P4,Point P5,Point P6,Point P7,Point P8,Point P9)
		{ n=11;p=new Point[n];p[0]=P0;p[1]=P1;p[2]=P2;p[3]=P3;p[4]=P4;p[5]=P5;p[6]=P6;p[7]=P7;p[8]=P8;p[9]=P9;p[10]=P0;}
		public Polygon(Point[] a,int 頂点数){ n=頂点数+1;if(n>0) p=new Point[n];for(int i=0;i<n;i++)p[i]=a[i];p[n-1]=a[0];}
		public Polygon(Line a)
		{ 
			if(a.p[0]==a.p[a.n-1])
			{ n=a.n;if(n>0){ p=new Point[n];for(int i=0;i<n;i++)p[i]=a.p[i];}}
			else
			{   
				n=a.n+1;
				if(n>0)	{ p=new Point[n];for(int i=0;i<n-1;i++)p[i]=a.p[i];	p[n-1]=a.p[0];}
			}
		}
		public Polygon(Polygon b)
		{
			n=b.n;if(n>0) p=new Point[n];for(int i=0;i<n;i++)p[i]=b.p[i];
		}

	}


	public class TMatrix : Point
	{
		public float[,] t=new float[4,4];
		public TMatrix(){ init();}
		public TMatrix
			(float D11,float D12,float D13,
			 float D21,float D22,float D23,
			 float D31,float D32,float D33)
		{
			t[0,0]=D11;t[0,1]=D12;t[0,2]=D13;
			t[1,0]=D21;t[1,1]=D22;t[1,2]=D23;
			t[2,0]=D31;t[2,1]=D32;t[2,2]=D33;
		}
		public TMatrix
			(float D11,float D12,float D13,float D14,
			 float D21,float D22,float D23,float D24,
			 float D31,float D32,float D33,float D34,
			 float D41,float D42,float D43,float D44)
		{
			t[0,0]=D11;t[0,1]=D12;t[0,2]=D13;t[0,3]=D14;
			t[1,0]=D21;t[1,1]=D22;t[1,2]=D23;t[1,3]=D24;
			t[2,0]=D31;t[2,1]=D32;t[2,2]=D33;t[2,3]=D34;
			t[3,0]=D41;t[3,1]=D42;t[3,2]=D43;t[3,3]=D44;
		}
		void init()
		{
			for(int i=0;i<4;i++)for(int j=0;j<4;j++)t[i,j]=0;
			for(int i=0;i<4;i++)t[i,i]=1;
		}
		public static Point operator *(TMatrix a,Point b)
		{
			return new Point(
				a.t[0,0]*b.X+a.t[0,1]*b.Y+a.t[0,2]*b.Z+a.t[0,3],
				a.t[1,0]*b.X+a.t[1,1]*b.Y+a.t[1,2]*b.Z+a.t[1,3],
				a.t[2,0]*b.X+a.t[2,1]*b.Y+a.t[2,2]*b.Z+a.t[2,3]);
		}
		public static Line operator *(TMatrix a,Line b)
		{
			Line c=new Line(b.n);
			for(int i=0;i<c.n;i++) c.p[i]=a*b.p[i];
			return c;
		}
		public static Polygon operator *(TMatrix a,Polygon b)
		{
			Polygon c=new Polygon(b.n-1);
			for(int i=0;i<c.n;i++) c.p[i]=a*b.p[i];
			return c;
		}
		public static TMatrix operator *(TMatrix a,TMatrix b)
		{
			TMatrix c=new TMatrix();
			for(int i=0;i<4;i++)
				for(int j=0 ;j<4;j++)
				{
					c.t[i,j]=0;
					for(int k=0;k<4;k++) c.t[i,j] += a.t[i,k]*b.t[k,j];
				}
			return c;

		}

		public static TMatrix scale(float scale_x, float scale_y, float scale_z)
		{
			TMatrix m=new TMatrix();
			m.t[0,0]=scale_x; m.t[1,1]=scale_y; m.t[2,2]=scale_z; return m;
		}
		public static TMatrix scale(Point P){ return scale(P.X, P.Y, P.Z);}
		public static TMatrix scale(float S){ return scale(S, S, S);}
		public static TMatrix operator *(TMatrix a,float S){return a*scale(S);}
		public static TMatrix operator *(float S,TMatrix a){return scale(S)*a;}
		public TMatrix perspect(float d)
		{
			TMatrix m=new TMatrix();
			m.t[2,3]=-d; m.t[3,2]=-1/d; return m;
		}
	}
	public class Surface  
	{
		public int n;
		public Polygon[] p;
		public Surface(int a){ n=a; if(n>0) p=new Polygon[n];}
		public Surface(Polygon a) {n=1; p=new Polygon[n];p[0]=a;}
		public Surface(Surface a) {n=a.n;if (n>0) p=new Polygon[n];for(int i=0;i<n; i++)p[i]=a.p[i];}
		public static Surface operator *(TMatrix a,Surface b)
		{ Surface c=new Surface(b.n); for(int i=0; i<c.n;i++) c.p[i]=a*b.p[i];return c;}
	}
	public class Ray
	{
		public Point o,d;
		public Ray (Point O, Point D){o=O;d=D;}
	}
	public class Ball:Point
	{
		public Point o;
		public float r;
		public Ball(Point O, float R){o=O;r=R;}
		public Ball(Point O){o=O;r=0;}
		public Ball(){o=new Point(0,0,0);r=0;}
		public float hit(Ray VR,ref Point P,ref Point N)
		{
			double ds=VR.d*(VR.o-o);
			double d2=ds*ds-(VR.o-o)*(VR.o-o)+r*r;
			if(d2<0) return 1E+8F;
			double t= -ds - Math.Sqrt(d2);
			//if(t<0.5) t= -ds + Math.Sqrt(d2);
			P=VR.o+(float)t * VR.d; N=unit(P-o);
			return (float)t;
		}
		public float hitV(Ray VR,ref Point P,ref Point N)
		{
			 double ds=VR.d*(VR.o-o);
			 double d2=ds*ds-(VR.o-o)*(VR.o-o)+r*r;
			 if(d2<0) return 1E+8F;
			 return (float)( -ds - Math.Sqrt(d2));
		 }

	}
	public class Wall:Point
	{
        public int n;
		public Polygon p;
		public Point nv;
		
		public Wall(Polygon pi)
		{p=new Polygon(pi); n=p.n;nv=unit((p.p[0]-p.p[2])%(p.p[1]-p.p[3]));}
		public float hit(Ray VR,ref Point P,ref Point N)
		{
			float t=nv*(p.p[1]-VR.o)/(nv*VR.d);  //距離
			if(t<=0) return 1E+8F;
			P=VR.o+t*VR.d; N=nv;
			for(int i=0; i<n-1; i++)
				if(N*((p.p[i+1]-p.p[i])%(P-p.p[i]))<0)return 1E+8F;
			return t;
		}

	}
	public class Attri:Point
	{
		public Color c;
		public Color[] tx;                 // 色，模様
		public float kd,ks,ke;             // 拡散反射率, 鏡面反射率, 環境光係数
		public float kr,kn,kt;             // 反射率, 屈折率, 透過率
		public int wi, hi;                 // ウインドウサイズ
		public int[] rn;
		public Attri(Color C, float D, float S, float E, float R, float N, float T)
		{  c=C; kd=D; ks=S; ke=E; kr=R; kn=N; kt=T; wi=0; tx=new Color[0]; rn=new int[1];rn[0]=1;}
		public Attri(Color C, float D, float S, float E, float R, float N)
		{  c=C; kd=D; ks=S; ke=E; kr=R; kn=N; kt=1; wi=0; tx=new Color[0]; rn=new int[1];rn[0]=1;}
		public Attri(Color C, float D, float S, float E, float R)
		{  c=C; kd=D; ks=S; ke=E; kr=R; kn=0; kt=1; wi=0; tx=new Color[0]; rn=new int[1];rn[0]=1;}
		public Attri(Color C, float D, float S, float E)
		{  c=C; kd=D; ks=S; ke=E; kr=0; kn=0; kt=1; wi=0; tx=new Color[0]; rn=new int[1];rn[0]=1;}
		public Attri(Color C, float D, float S)
		{  c=C; kd=D; ks=S; ke=0.2F; kr=0; kn=0; kt=1; wi=0; tx=new Color[0]; rn=new int[1];rn[0]=1;}
		public Attri(Color C, float D)
		{  c=C; kd=D; ks=0.8F; ke=0.2F; kr=0; kn=0; kt=1; wi=0; tx=new Color[0]; rn=new int[1];rn[0]=1;}
		public Attri(Color C)
		{  c=C; kd=0.7F; ks=0.8F; ke=0.2F; kr=0; kn=0; kt=1; wi=0; tx=new Color[0]; rn=new int[1];rn[0]=1;}
		public Attri()
		{  c=Color.White; kd=0.7F; ks=0.8F; ke=0.2F; kr=0; kn=0; kt=1; wi=0; tx=new Color[0]; rn=new int[1];rn[0]=1;}
		public Attri(Attri a)
		{  c=a.c; kd=a.kd; ks=a.ks; ke=a.ke; kr=a.kr; kn=a.kn; kt=a.kt; wi=a.wi; hi=a.hi; tx=a.tx; rn=a.rn;}
		public Color shading(Point V, Point L, Point N, Attri A, double s)
		{
			Point RV=L-2*(L*N)*N;			
			double D1=(-N*L)/(len(N)*len(L));　 // 拡散および環境光
			if(D1<0)D1=0;
			D1=A.kd*s*D1+A.ke;　                         
			
			double D2=s*A.ks*((-RV*V))/(len(RV)*len(V));　　　   // 反射
			if(D2<0)D2=0;
			int R=(int)(D1*(double)(A.c.R)+D2*255);if(R>255)R=255;
			int G=(int)(D1*(double)(A.c.G)+D2*255);if(G>255)G=255;
			int B=(int)(D1*(double)(A.c.B)+D2*255);if(B>255)B=255;
			return Color.FromArgb(R,G,B);
		}

	}

	public class Quad:Point    //二次局面
	{
		public int type;             //次数
		public TMatrix A;            //2次係数
		public Point   B;            //1次係数
		public float   C;            //0次係数
		public Quad() { type=0; C=0;A=new TMatrix();B=new Point(0,0,0); }
		public Quad(Point s2, Point s1, float s0)
		{
			A=scale(s2);B=s1;C=s0;
			if     (!(s2==new Point(0,0,0))) type=2;
			else if(!(s1==new Point(0,0,0))) type=1;
			else                           type=0;
		}
		public bool inner(Point P)
		{
			if(type==1)if((P*B + C)<=0) return true;else return false;
			else if((A*P)*P+this.B*P+C<=0) return true;
			return false;
		}
		public float hit(Ray VR, ref Point P, ref Point N)
		{
			Point V0 = VR.o;
			Point V  = VR.d;
			float t  = 1E8F; double ad = 0; 
			if(type ==0) return t;
			if(type ==1) { t=-(B*V0+C)/(B*V);}
			else
			{
				double a = (A*V)*V, b=(2*A*V0+B)*V, c=(A*V0+B)*V0+C;
				double d = b*b-4*a*c;
				if(d<0) return 1E8F;
				ad=Math.Sqrt(d)/Math.Abs(2*a);
				t=(float)(-b/(2*a)-ad);
				if(t<0.5F) t=t+2*(float)ad;
			}
			if(t<0.5F) return 1E8F;
			P=VR.o+t*VR.d;
			N=unit(2*A*P+B);
			return t;

		}
		public static Quad operator *(TMatrix a,Quad b)   
		{
			Quad c=b;
			TMatrix r=inverse(a);
			Point o=-(a*(new Point(0,0,0)));
			r.t[0,3]=r.t[1,3]=r.t[2,3]=0;
			c.A=trans(r)*b.A*r; c.B=trans(r)*b.B+2*c.A*o;
			c.C=b.C+(c.B-c.A*o)*o;
			return c;
		}
		public static Quad operator -(Quad b)   
		{
			Quad c=b; c.A=-1*b.A;c.B=-b.B; c.C=-b.C; return c;
		}
		static TMatrix scale(float scale_x, float scale_y, float scale_z)
		{
			TMatrix m=new TMatrix();
			m.t[0,0]=scale_x; m.t[1,1]=scale_y; m.t[2,2]=scale_z; return m;
		}
		static TMatrix scale(Point P){ return scale(P.X, P.Y, P.Z);}
		static TMatrix scale(float S){ return scale(S, S, S);}
		static TMatrix inverse(TMatrix a)
		{
			float d=0;
			for(int i=0; i<3; i++)
				d = d + a.t[0,i]*a.t[1, (i+1)%3]*a.t[2, (i+2) % 3]
					- a.t[0,i]*a.t[1, (i+2)%3]*a.t[2, (i+1) % 3];
			TMatrix b=new TMatrix();
			if(d != 0)
				for(int i=0; i<3; i++)
					for(int j=0; j<3; j++)
						b.t[i,j]=(a.t[(j+1)%3,(i+1)%3]*a.t[(j+2)%3, (i+2)%3]
							-a.t[(j+1)%3,(i+2)%3]*a.t[(j+2)%3, (i+1)%3])/d;
			for(int i=0; i<3; i++) b.t[i,3]=-a.t[i,3];
			return b;

		}
		static TMatrix trans(TMatrix a)
		{
			TMatrix c=new TMatrix();
			for(int i=0; i<4; i++)
				for(int j=0; j<4; j++) c.t[i,j]=a.t[j,i];
			return c;
		}
	}
	public class Csgm :Point  //Constructive Solid Geometry Model
	{
		private const int   AND     = 1;
		private const int   OR      = 2;
		private const int   MINUS   = 3;
		private const float 交点なし= 1E8F;

		public Csgm c1,c2;
		public Quad h1;
		public int op;
		public Ball BV;
		public Csgm(){ op=0;h1=new Quad();}
		public Csgm(Point s2,Point s1,float s0)	{ h1=new Quad(s2,s1,s0);op=0;}
		public Csgm(Csgm a,Csgm b, int md)
		{
			BV=new Ball(new Point(0,0,0),-1);
			c1=new Csgm(); c2=new Csgm();
			if(md==MINUS){c2 =-b; op=AND;}
			else         {c2 = b;  op=md;}
			c1=a;
		}
		public Csgm(Csgm a,Csgm b){ new Csgm(a,b,0);}
		public Csgm(Csgm a)
		{
			BV=a.BV;
			c1=new Csgm();c2=new Csgm();h1=new Quad();
			if(a.op==0)
			{
				c1=a.c1;c2=a.c2;
			}
			else h1=a.h1;
			op=a.op;
		}

		public Csgm(Csgm a, Ball bv)
		{
			BV=bv;
			c1=new Csgm();c2=new Csgm();h1=new Quad();
			if(op==0)
			{
				c1=a.c1;c2=a.c2;
			}
			else h1=a.h1;
			op=a.op;
		}

		public static Csgm operator *(TMatrix a,Csgm b)
		{
			Csgm  c=b;
			if(b.op==0){ c.h1=a*(b.h1);return c;}
			c.c1=a*(b.c1);c.c2=a*(b.c2); return c;
		}
		public static Csgm operator -(Csgm b)
		{
			Csgm  c=b;
			if(b.op==0){ c.h1=-(b.h1);return c;}
			if(b.op==AND) c.op=OR;
			else if(b.op==OR) c.op=AND;
            c.c1=-(b.c1);c.c2=-(b.c2); return c;
		}
		public static Csgm operator +(Csgm a, Csgm b){return new Csgm(a,b,OR);}
		public static Csgm operator -(Csgm a, Csgm b){return new Csgm(a,b,MINUS);}
		public static Csgm operator *(Csgm a, Csgm b){return new Csgm(a,b,AND);}
		public bool inner(Point P)
		{
			if(op==AND) return c1.inner(P) & c2.inner(P);
			if(op==OR)  return c1.inner(P) | c2.inner(P);
			return h1.inner(P);
		}
		public bool outer(Point P){ return !inner(P);}
		public float hit(Ray VR,ref Point P, ref Point N)
		{
			if((BV.r>0) && (BV.hit(VR,ref P,ref N)>=交点なし)) return 交点なし;
			Point P1=new Point(),P2=new Point(),N1=new Point(),N2=new Point();
			double t1,t2;
			if(op==AND)
			{
				Ray V=VR;
				for(int i=0;i<10;i++)
				{
					t1=c1.hit(V, ref P1,ref N1);
					t2=c2.hit(V, ref P2,ref N2);
					if(t1>=交点なし && t2>=交点なし) return 交点なし;
					if(Math.Abs(t1-t2)<0.5)
					{
						Point PP=V.o+(float)(Math.Max(t1,t2)+0.5)*V.d;
						if(c1.inner(PP) && c2.inner(PP))
						{ P=V.o+(float)(Math.Max(t1,t2))*V.d; N=N1; return len(VR.o-P);}
						else return 交点なし;
					}
					if(t1<t2)
						if(c2.inner(P1)) {N=N1;P=P1; return len(VR.o-P1);}
						else             {V.o=P1;}
					else if(c1.inner(P2)){N=N2;P=P2; return len(VR.o-P2);}
					else                 {V.o=P2;}
				}
				return 交点なし;
			}
			else if(op==OR)
			{
				Ray V=VR;
				for(int i=0;i<10;i++)
				{
					t1=c1.hit(V, ref P1,ref N1);
					t2=c2.hit(V, ref P2,ref N2);
					if(t1>=交点なし && t2>=交点なし) return 交点なし;
					if(t1<t2)
						if(c2.outer(P1)) {N=N1;P=P1; return len(VR.o-P1);}
						else             {V.o=P1;}
					else if(c1.outer(P2)){N=N2;P=P2; return len(VR.o-P2);}
					else                 {V.o=P2;}
				}
			}
			else if(op==0) return h1.hit(VR,ref P,ref N);
			return 交点なし;
		}
		public Csgm ball(Point s)
		{
			Point s2=new Point((s.X==0)?0:1/(s.X*s.X),(s.Y==0)?0:1/(s.Y*s.Y),(s.Z==0)?0:1/(s.Z*s.Z));
			return new Csgm(s2, new Point(0,0,0),-1);
		}
		public Csgm ball    (float r)          { return ball(new Point(r,r,r));}
		public Csgm cylinder(Point s)          { return ball(s);}
		public Csgm plane   (Point a, float b) { return new Csgm(new Point(0, 0, 0), a, b);}
		public Csgm plane   (Point a)          { return plane(a,0);}
		public Csgm plane_xy()                 { return new Csgm(new Point(0, 0, 0), new Point(0, 0, 1),0);}
		public Csgm plane_yz()                 { return new Csgm(new Point(0, 0, 0), new Point(1, 0, 0),0);}
		public Csgm plane_xz()                 { return new Csgm(new Point(0, 0, 0), new Point(0, 1, 0),0);}
		public Csgm boad_xy(float b)
		{  
			return (new Csgm(new Point(0, 0, 0), new Point(0, 0, -1),-b/2))*
			       (new Csgm(new Point(0, 0, 0), new Point(0, 0,  1),-b/2));
		}
	}
	public class Obj
	{
		private const int WALL    =0x10;
		private const int BALL    =0x20;
		private const int CSGM    =0x30;
		private const int FREE    =0x40;
		private const int BVOL    =0xF0;

		public Wall  w;
		public Ball  b;
		public Csgm  c;
		public int   type;
		public Attri att;
		public int   ele;

		public Obj(Wall W, Attri a){w=W;type=WALL; att=a;} 
		public Obj(Ball B, Attri a){b=B;type=BALL; att=a;} 
		public Obj(Csgm C, Attri a){c=C;type=CSGM; att=a;} 
	}
	public class Scene: Attri
	{
		private const int MAX_OBJ =1024;
		private const int WALL    =0x10;
		private const int BALL    =0x20;
		private const int CSGM    =0x30;
		private const int FREE    =0x40;
		private const int BVOL    =0xF0;
		private const float 交点なし= 1E8F;

		public Obj   [] obj  = new Obj   [MAX_OBJ];//物体
		public int    n;                           //物体の個数
		public Color  bgc;                         //背景色
		public bool   md;                          //表示の有無
		public Scene(Color C, bool M){ n=0; md = M   ; bgc = C;}
		public Scene(Color C)        { n=0; md = true; bgc = C;}
		public Scene()               { n=0; md = true; bgc = Color.Black;}
		public Color bgcolor(Ray VR) {return bgc;}
		public void regist(Wall a,Attri c){ obj[n++]=new Obj(a,c); 	}
		public void regist(Ball a,Attri c){ obj[n++]=new Obj(a,c);	}
		public void regist(Csgm a,Attri c, Ball bv)
		{ Csgm X =new Csgm(a,bv); obj[n++]=new Obj(X,c);	}
		public void regist(Csgm a,Attri c)
		{ Csgm X =new Csgm(a,new Ball(new Point(0,0,0),-1)); obj[n++]=new Obj(X,c);	}
		public Attri attri(int i,Point P,Point N)
		{   
			//MessageBox.Show("obj["+i+"]=" + obj[i].att.ToString());
			return obj[i].att;
		}
		public int hit(Ray VR, ref Point P, ref Point N, ref int hi)
		{
			float t=交点なし, tt=0F;
			Point HP=new Point(0,0,0),HN=new Point(0,0,0);
			hi=-1;
			for(int i=0;i<n; i++)
			{
				if     (obj[i].type==WALL)tt=obj[i].w.hit(VR,ref HP,ref HN);
				else if(obj[i].type==BALL)tt=obj[i].b.hit(VR,ref HP,ref HN);
				else if(obj[i].type==BALL)tt=obj[i].b.hit(VR,ref HP,ref HN);
				if(tt>0.5 && tt<t){t=tt; hi=i; P=HP; N=HN;}
			}
			return hi;
		}
		public float trans(Ray VR,float L) //照明光の強さ
		{
			float s=1F, tt=0F;
			Point HP=new Point(0,0,0),HN=new Point(0,0,0);
			for(int i=0;i<n; i++)
			{
				if     (obj[i].type==WALL)tt=obj[i].w.hit(VR,ref HP,ref HN);
				else if(obj[i].type==BALL)tt=obj[i].b.hit(VR,ref HP,ref HN);
				else if(obj[i].type==BALL)tt=obj[i].b.hit(VR,ref HP,ref HN);
				if(tt>0.5 && tt<L-0.5)
					if(obj[i].att.kn<=0)s=0; else s= s*(1-obj[i].att.kr);
				if(s<=0) break;
			}
			return s;
		}
		public Color color(Ray VR, Point L0,int rn)
		{
			int r=rn; return color(VR,ref L0,ref r,1);

		}
		public Color color(Ray VR, ref Point L0,ref int rn, float am)
		{
			if((rn==0) || (am<0.01)) return Color.Black;
			Point Ps=new Point(0,0,0), N=new Point(0,0,0);
			int hi=-1, rnt=rn-1,rnr=rn-1;   // hi: 視線が交わった物体の番号
			if(hit(VR,ref Ps,ref N,ref hi)<0){ rn=0; return bgcolor(VR);}

			Color Cs=Color.Black,Cr=Color.Black,Ct=Color.Black;
			//MessageBox.Show("hi=" +hi.ToString());
			Attri A  = attri(hi, Ps,N);
			float kr = A.kr;
			Point L  = unit(Ps-L0);                //照明光方向ベクトル
			Ray   LR = new Ray(Ps,-L);             //照明光線定義
			float s  = trans(LR, len(L0-Ps));      //交点が日陰
			Cs=shading(VR.d,L,N,A,s);              //交点自身の色
			if(A.kn>0)
			{                                      //屈折光処理
				float n=A.kn;                      
				if(N*VR.d>0) {n=1/n; N=-N;}        //物体から出る場合
				float c = -VR.d*N, g2 = n*n+c*c-1;
				if(g2>0)                           //全反射でない場合
				{
					float g=(float)(Math.Sqrt((double)g2));
					Ray   TR = new Ray(Ps, (VR.d+(c-g)*N)/n);
					double cmg=c-g, cpg=c+g, nmg=n*n*c-g, npg=n*n*c+g;
					kr = (float)((cmg*cmg/(cpg*cpg)+nmg*nmg/(npg*npg))/2);
					Ct = color(TR,ref L0,ref rnt, am*(1-kr)*A.kt);
				}
				else kr=1;                         //全反射
    		}
			if(kr>0)                               
			{
				if(N*VR.d>0) N=-N;
				Ray RR=new Ray(Ps,VR.d-2*(VR.d*N)*N); // 反射視線定義
				Cr=color(RR,ref L0,ref rnr,am*kr);    // 反射視点から見える物体色
			}
			rn=((rnr>rnt)?rnr:rnt)+1;
			return mixColor(Cr,Ct,Cs,kr,(1-kr)*A.kt,1);
    	}
		private Color mixColor(Color c1,Color c2,Color c3,float a1,float a2,float a3)
		{
			float R1=(float)c1.R, G1=(float)c1.G, B1=(float)c1.B;
			float R2=(float)c2.R, G2=(float)c2.G, B2=(float)c2.B;
			float R3=(float)c3.R, G3=(float)c3.G, B3=(float)c3.B;
			int R=(int)(R1*a1+R2*a2+R3*a3); if(R<0)R=0;else if(R>255) R=255;
			int G=(int)(G1*a1+G2*a2+G3*a3); if(G<0)G=0;else if(G>255) G=255;
			int B=(int)(B1*a1+B2*a2+B3*a3); if(B<0)B=0;else if(B>255) B=255;
			return Color.FromArgb(R,G,B);
		}	
	}
	public class Hyper:Point
	{
		private const float 交点なし= 1E8F;

		public Point S1,S2; public float S0;
		public Hyper(Point s2, Point s1, float s0){ S2=s2; S1=s1;S0=s0;}
		public float hit(Ray VR, ref Point P,ref Point N)
		{
			Point V0 = VR.o, V  = VR.d;
			TMatrix scaleS2=scale(S2);
			Point scaleS2_V=scaleS2*V;
			float d=(scaleS2_V*V0+0.5F*S1*V)*(scaleS2_V*V0+0.5F*S1*V)-
				    (scaleS2_V*V)*((scaleS2*V0)*V0+S1*V0+S0);
			if(d<0) return 交点なし;
			float ad=(float)(Math.Sqrt((double)d)/Math.Abs((double)(scaleS2_V*V)));
			float t=-(scaleS2_V*V0+0.5F*S1*V)/(scaleS2_V*V)-ad;
			P=VR.o+t*VR.d;
			if(t<0||P.X<-120F||120F<P.X||P.Y<-120F||120F<P.Y||P.Z<-120F||120F<P.Z) t=t+2*ad;
			P=VR.o+t*VR.d;
			if(t<0||P.X<-120F||120F<P.X||P.Y<-120F||120F<P.Y||P.Z<-120F||120F<P.Z) return 交点なし;
			N=unit(2*scale(S2)*P+S1);
			return t;
		}
		public static TMatrix scale(float scale_x, float scale_y, float scale_z)
		{
			TMatrix m=new TMatrix();
			m.t[0,0]=scale_x; m.t[1,1]=scale_y; m.t[2,2]=scale_z; return m;
		}
		public static TMatrix scale(Point P){ return scale(P.X, P.Y, P.Z);}
		public static TMatrix scale(float S){ return scale(S, S, S);}

	}
}