(* ::Package:: *)

(* DyGeom.m   by Kazushi Ahara *)
(* pre-alpha version 0.03 *)
If[KETPIC5, <<"ketpic5.m"; ,<<"ketpic6.m"; ](* You need ketpic6.m for mathematics *)
NextPointName = "P";
NextPointNumber = 0;
NextLineName = "L";
NextLineNumber = 0;
NextCircleName = "C";
NextCircleNumber = 0;
NextMarkName = "M";
NextMarkNumber = 0;
MyTeXFileName="zu.tex";

GetNewPointName[] := Module[{},
	NextPointNumber=NextPointNumber+1;
	NextPointName<>ToString[NextPointNumber]
]

GetNewLineName[] := Module[{},
	NextLineNumber ++;
	NextLineName<>ToString[NextLineNumber]
]

GetNewCircleName[] := Module[{},
	NextCircleNumber ++;
	NextCircleName<>ToString[NextCircleNumber]
]

GetNewMarkName[] := Module[{},
	NextMarkNumber ++;
	NextMarkName<>ToString[NextMarkNumber]
]

ReduceTriple[arg_]:=Module[{myGCD,ret},
If[IntegerQ[arg[[1]]] && IntegerQ[arg[[2]]]&&IntegerQ[arg[[3]]], 
	myGCD=GCD[arg[[1]],arg[[2]],arg[[3]]];
	ret=If[myGCD!=0, arg/myGCD, arg]; , 
	ret=arg;
];
ret
]

ReduceTriplePolynomial[arg_]:=Module[{myGCD,newLL},
newLL=arg;
myGCD=PolynomialGCD[arg[[1]],arg[[2]],arg[[3]]];
newLL=Together[newLL/myGCD];
newLL
]

FreePoint[argseq__] := Module[ {args,ptname,ptx,pty},
args={argseq};
ptname = GetNewPointName[];
If[NumberQ[args[[1]]], ptx=args[[1]], ptx = NextPointNumber];
args = Drop[args , 1];
If[NumberQ[args[[1]]], pty=args[[1]], pty = 0];
args = Drop[args , 1];
Append[{"Point",ptname,{ptx,pty,1},{ptname<>"x",ptname<>"y",ptname<>"z"}},args]
]

ConnectPointsCoord[argseg__] := Module[ {args},
args={argseg};
{ args[[1]][[2]]*args[[2]][[3]]-args[[1]][[3]]*args[[2]][[2]],
args[[1]][[3]]*args[[2]][[1]]-args[[1]][[1]]*args[[2]][[3]],
args[[1]][[1]]*args[[2]][[2]]-args[[1]][[2]]*args[[2]][[1]] }
]

ConnectPoints[argseq__] := Module[ {args,lnname,p1,p2,ret},
args={argseq};
lnname = GetNewLineName[];
ret={"Line",lnname,{0,0,0},{0,0,0},{}};
If[Length[args]>=2,
	p1=args[[1]];
	p2=args[[2]];
	args=Drop[args,2];
	If[p1[[1]]=="Point" && p2[[1]]=="Point",
		ret={"Line",lnname, 
		ReduceTriple[ConnectPointsCoord[p1[[3]],p2[[3]]]],
		ReduceTriplePolynomial[ConnectPointsCoord[p1[[4]],p2[[4]]]],
		args} ;
	];
];
ret
]

MidPointCoord[argseq__] := Module[ {args,p1,p2},
args={argseq};
p1=args[[1]];
p2=args[[2]];
{p1[[1]]*p2[[3]]+p2[[1]]*p1[[3]],
p1[[2]]*p2[[3]]+p2[[2]]*p1[[3]],
2*p1[[3]]*p2[[3]]}
]

MidPoint[argseq__] := Module[ {args,ptname,p1,p2,ret},
args={argseq};
ptname = GetNewPointName[];
ret={"Point",ptname,{0,0,0},{0,0,0},{}};
If[Length[args]>=2,
	p1=args[[1]];
	p2=args[[2]];
	args=Drop[args,2];
	If[p1[[1]]=="Point" && p2[[1]]=="Point",
		ret={"Point",ptname,
		ReduceTriple[MidPointCoord[p1[[3]],p2[[3]]]],
		MidPointCoord[p1[[4]],p2[[4]]],
		args}; 
	];
];
ret
]

IntersectionLL[argseq__]:=Module[{args,ptname,l1,l2,ret,myGCD,newLL},
args={argseq};
ptname=GetNewPointName[];
ret={"Point",ptname,{0,0,0},{0,0,0},{}};
	l1=args[[1]];
	l2=args[[2]];
	args=Drop[args,2];
	If[l1[[1]]=="Line" && l2[[1]]=="Line",
		ret={"Point",ptname, 
		ReduceTriple[ConnectPointsCoord[l1[[3]],l2[[3]]]],
		ReduceTriplePolynomial[ConnectPointsCoord[l1[[4]],l2[[4]]]],
		args} ;
	];
ret
]

ParallelLineCoord[argseq__] := Module[ {args,l1,p1},
args={argseq};
l1=args[[1]];
p1=args[[2]];
{l1[[1]]*p1[[3]],l1[[2]]*p1[[3]],-(l1[[1]]*p1[[1]]+l1[[2]]*p1[[2]])}
]

ParallelLine[argseq__] := Module[ {args,lnname,l1,p1,ret},
args={argseq};
lnname = GetNewLineName[];
ret={"Line",lnname,{0,0,0},{0,0,0},{}};
If[Length[args]>=2,
	l1=args[[1]];
	If[l1[[1]]=="Point", p1=l1;l1=args[[2]], p1=args[[2]];];
	args=Drop[args,2];
	If[l1[[1]]=="Line" && p1[[1]]=="Point",
		ret={"Line",lnname,
		ReduceTriple[ParallelLineCoord[l1[[3]],p1[[3]]]],
		ParallelLineCoord[l1[[4]],p1[[4]]],
		args}; 
	];
];
ret
]

PerpendicularCoord[argseq__] := Module[ {args,l1,p1},
args={argseq};
l1=args[[1]];
p1=args[[2]];
{-l1[[2]]*p1[[3]],l1[[1]]*p1[[3]],l1[[2]]*p1[[1]]-l1[[1]]*p1[[2]]}
]

PerpendicularLine[argseq__] := Module[ {args,lnname,l1,p1,ret},
args={argseq};
lnname = GetNewLineName[];
ret={"Line",lnname,{0,0,0},{0,0,0},{}};
If[Length[args]>=2,
	l1=args[[1]];
	If[l1[[1]]=="Point", p1=l1;l1=args[[2]], p1=args[[2]];];
	args=Drop[args,2];
	If[l1[[1]]=="Line" && p1[[1]]=="Point",
		ret={"Line",lnname,
		ReduceTriple[PerpendicularCoord[l1[[3]],p1[[3]]]],
		PerpendicularCoord[l1[[4]],p1[[4]]],
		args}; 
	];
];
ret
]

AngleLineCoord[argseq__] := Module[ {args,l1,p1,theta},
args={argseq};
l1=args[[1]];
p1=args[[2]];
theta=args[[3]];
{(Cos[theta]*l1[[1]]-Sin[theta]*l1[[2]])*p1[[3]],(Sin[theta]*l1[[1]]+Cos[theta]*l1[[2]])*p1[[3]],-((Cos[theta]*l1[[1]]-Sin[theta]*l1[[2]])*p1[[1]]+(Sin[theta]*l1[[1]]+Cos[theta]*l1[[2]])*p1[[2]])}
]

AngleLine[argseq__] := Module[ {args,lnname,l1,p1,theta,ret},
args={argseq};
lnname = GetNewLineName[];
ret={"Line",lnname,{0,0,0},{0,0,0},{}};
If[Length[args]>=3 && NumberQ[args[[3]]],
	l1=args[[1]];
	If[l1[[1]]=="Point", p1=l1;l1=args[[2]], p1=args[[2]];];
	theta=args[[3]];
	args=Drop[args,3];
	If[IntegerQ[theta], theta=theta*Pi/180];
	If[l1[[1]]=="Line" && p1[[1]]=="Point",
		ret={"Line",lnname,
		ReduceTriple[AngleLineCoord[l1[[3]],p1[[3]],theta]],
		AngleLineCoord[l1[[4]],p1[[4]],theta],
		args}; 
	];
];
ret
]

CircleCoord[argseq__] := Module[ {args,p1,p2},
args={argseq};
p1=args[[1]];
p2=args[[2]];
{p1[[3]]*p2[[3]]^2,
-2*p1[[1]]*p2[[3]]^2,
-2*p1[[2]]*p2[[3]]^2,
-p1[[3]]*(p2[[1]]^2+p2[[2]]^2)+2*(p1[[1]]*p2[[1]]+p1[[2]]*p2[[2]])*p2[[3]]}
]

CircleByCenter[argseq__]:=Module[{args,cnname,p1,p2,ret},
args={argseq};
cnname=GetNewCircleName[];
ret={"Circle",cnname,{0,0,0,0},{0,0,0,0},{}};
If[Length[args]>=2,
	p1=args[[1]];
	p2=args[[2]];
	args=Drop[args,2];
	If[p1[[1]]=="Point" && p2[[1]]=="Point",
		ret={"Circle",cnname,
		CircleCoord[p1[[3]],p2[[3]]],
		CircleCoord[p1[[4]],p2[[4]]],
		args}; 
	];
];
ret
]

BiangularBy3Coor[arg1_,arg2_,arg3_]:=Module[{d12,d23,pn},
d12=Sqrt[(arg1[[1]]/arg1[[3]]-arg2[[1]]/arg2[[3]])^2+(arg1[[2]]/arg1[[3]]-arg2[[2]]/arg2[[3]])^2];
d23=Sqrt[(arg3[[1]]/arg3[[3]]-arg2[[1]]/arg2[[3]])^2+(arg3[[2]]/arg3[[3]]-arg2[[2]]/arg2[[3]])^2];
pn={arg1[[1]]*arg3[[3]]*d23+arg3[[1]]*arg1[[3]]*d12,arg1[[2]]*arg3[[3]]*d23+arg3[[2]]*arg1[[3]]*d12,(d12+d23)*arg1[[3]]*arg3[[3]]};
ConnectPointsCoord[pn,arg2]
]

BiangularBy3[argseq__]:=Module[{args,lnname,p1,p2,p3,pn,len12,len23,ret},
args={argseq};
lnname = GetNewLineName[];
ret={"Line",lnname,{0,0,0},{0,0,0},{}};
If[Length[args]>=3,
	p1=args[[1]];
	p2=args[[2]];
	p3=args[[3]];
	args=Drop[args,3];
	If[p1[[1]]=="Point" && p2[[1]]=="Point" && p3[[1]]=="Point" && p1[[3]][[3]]!=0 && p2[[3]][[3]]!=0 && p3[[3]][[3]]!=0,
		ret={"Line",lnname,
		N[BiangularBy3Coor[p1[[3]],p2[[3]],p3[[3]]]],
		BiangularBy3Coor[p1[[4]],p2[[4]],p3[[4]]],
		args};
	];
];
ret
]

AngleMarkBy3[argseq__]:=Module[{args,mkname,p1,p2,p3,ret},
args={argseq};
mkname = GetNewMarkName[];
ret={"Mark",mkname,{0,0,0},{0,0,0},{}};
If[Length[args]>=3,
	p1=args[[1]];
	p2=args[[2]];
	p3=args[[3]];
	args=Drop[args,3];
	If[p1[[1]]=="Point" && p2[[1]]=="Point" && p3[[1]]=="Point" && p1[[3]][[3]]!=0 && p2[[3]][[3]]!=0 && p3[[3]][[3]]!=0,
		ret={"Mark",mkname,
		{p1[[3]],p2[[3]],p3[[3]]},{0,0,0},
		args};
	];
];
ret
]

AngleMarkPosotion[arg1_,arg2_,arg3_,mkpos_]:=Module[{d12,d32,d42,arg4},
d12=Norm[arg1-arg2];
d32=Norm[arg3-arg2];
arg4=(d32*arg1+d12*arg3)/(d12+d32);
d42=Norm[arg4-arg2];
arg2+(arg4-arg2)*mkpos/d42
]

AddOption[arg1_,arg2_]:=Module[{arg15},
arg15=Append[arg1[[5]],arg2];
Append[Drop[arg1,-1],arg15]
]

GetXYdomain[argseq__]:= Module[ {obj,args,xmax,xmin,ymax,ymin,r2},
args={argseq};
ymax=xmax=-999;
ymin=xmin=999;
For[i=1,i<=Length[args],i++,
	obj=args[[i]][[3]];
	If[args[[i]][[1]]=="Point" && obj[[3]]!=0 , 
		If[xmax<obj[[1]]/obj[[3]],xmax=obj[[1]]/obj[[3]]];
		If[xmin>obj[[1]]/obj[[3]],xmin=obj[[1]]/obj[[3]]];
		If[ymax<obj[[2]]/obj[[3]],ymax=obj[[2]]/obj[[3]]];
		If[ymin>obj[[2]]/obj[[3]],ymin=obj[[2]]/obj[[3]]]
	];
	If[args[[i]][[1]]=="Circle" && obj[[1]]!=0 , 
		r2=-obj[[4]]/obj[[1]]+(obj[[2]]^2+obj[[3]]^2)/(4*obj[[1]]^2);
		If[r2>0 && xmax < (-obj[[2]]/(2*obj[[1]]))+Sqrt[r2], xmax = (-obj[[2]]/(2*obj[[1]]))+Sqrt[r2]];
		If[r2>0 && xmin > (-obj[[2]]/(2*obj[[1]]))-Sqrt[r2], xmin = (-obj[[2]]/(2*obj[[1]]))-Sqrt[r2]];
		If[r2>0 && ymax < (-obj[[3]]/(2*obj[[1]]))+Sqrt[r2], ymax = (-obj[[3]]/(2*obj[[1]]))+Sqrt[r2]];
		If[r2>0 && ymin > (-obj[[3]]/(2*obj[[1]]))-Sqrt[r2], ymin = (-obj[[3]]/(2*obj[[1]]))-Sqrt[r2]];
	];	
	If[xmax>=xmin+1, XMAX=(xmax*10-xmin)/9;XMIN=(xmin*10-xmax)/9;,
	XMAX=xmax+1;XMIN=xmin-1;];
	If[ymax>=ymin+1, YMAX=(ymax*10-ymin)/9;YMIN=(ymin*10-ymax)/9;,
	YMAX=ymax+1;YMIN=ymin-1;];
];
]

AlmostZero[arg_]:=Module[{},
If[arg==0, True, False]
]

Incidency[arg1_,arg2_]:=Module[{txt,testfunction},
txt="";
If[arg1[[1]]=="Point" && arg2[[1]]=="Line",
	If[Abs[arg1[[3]].arg2[[3]]]>0.0000001, txt="No incidency";,
		testfunction=Factor[arg1[[4]].arg2[[4]]];
		If[AlmostZero[testfunction], 
			txt=arg1[[2]]<>" is on the line "<>arg2[[2]]<>"."; , 
			txt="No incidency";
		];
	];
];
If[arg1[[1]]=="Line" && arg2[[1]]=="Point",
	If[Abs[arg1[[3]].arg2[[3]]]>0.0000001, txt="No incidency";,
		testfunction=Factor[arg1[[4]].arg2[[4]]];
		If[AlmostZero[testfunction], 
			txt=arg2[[2]]<>" is on the line "<>arg1[[2]]<>"."; , 
			txt="No incidency";
		];
	];
];
txt
]

MakePlotData[argseq__]:=Module[{args,pd,pointdata,textdata,obj,xsect1,xsect2,option,ptname,
ptposition,ptsize,radius,clip1,clip2,x1,x2,y1,y2,r2,xxmax,xxmin,yymax,yymin,d12,d23},
args={argseq};
pd={};
pointdata={};
textdata={};
For[i=1,i<=Length[args],i++,
	obj=args[[i]][[3]];
	If[args[[i]][[1]]=="Point" && obj[[3]]!=0 , 
		ptname="";ptposition="";ptsize="";
		option=args[[i]][[5]];
		(*Print[option];*)
		For[j=1,j<=Length[option],j++,
			Switch[option[[j]][[1]], 
				pointname, ptname=option[[j]][[2]],
				pointposition, ptposition=option[[j]][[2]],
				pointsize, ptsize=option[[j]][[2]]
			];
		];
		If[ptname!="",
			textdata=Append[textdata,{obj[[1]]/obj[[3]],obj[[2]]/obj[[3]],ptposition,ptname}]
		];
		If[NumberQ[ptsize],
			radius=ptsize, 
			radius=Switch[ptsize,"small",1,"midium",2,"large",3,_,2]
		];
		pointdata=Append[pointdata,{obj[[1]]/obj[[3]],obj[[2]]/obj[[3]],radius}];
	];
	If[args[[i]][[1]]=="Line", 
		option=args[[i]][[5]];
		clip1=clip2={};
		For[j=1,j<=Length[option],j++,
			Switch[option[[j]][[1]], 
				clipping, clip1=option[[j]][[2]][[1]];clip2=option[[j]][[2]][[2]]
			];
		];
		xxmax=XMAX;xxmin=XMIN;yymax=YMAX;yymin=YMIN;
		(*Print["Here I am. ",Length[clip1]," and ",Length[clip2]];*)
		If[Length[clip1]>=3 && Length[clip2]>=3, 
			(*Print["I'm in depth 2. ",clip1[[3]][[3]]];*)
			If[clip1[[1]]=="Point" && clip2[[1]]=="Point" && clip1[[3]][[3]]!=0 && clip2[[3]][[3]]!=0,
				Print[clip1];
				xxmax=yymax=-999; xxmin=yymin=999;
				x1=clip1[[3]][[1]]/clip1[[3]][[3]];
				If[xxmax<x1,xxmax=x1];
				If[xxmin>x1,xxmin=x1];
				x2=clip2[[3]][[1]]/clip2[[3]][[3]];
				If[xxmax<x2,xxmax=x2];
				If[xxmin>x2,xxmin=x2];
				y1=clip1[[3]][[2]]/clip1[[3]][[3]];
				If[yymax<y1,yymax=y1];
				If[yymin>y1,yymin=y1];
				y2=clip2[[3]][[2]]/clip2[[3]][[3]];
				If[yymax<y2,yymax=y2];
				If[yymin>y2,yymin=y2];
			];
			(*Print[xxmax,xxmin,yymax,yymin]; *)
		];
		If[obj[[2]]==0 && obj[[1]]!=0,
			y1=yymax;
			y2=yymin;
			pd=Append[pd,listplotdata[{{-obj[[3]]/obj[[1]],y1},{-obj[[3]]/obj[[1]],y2}}]]
		];
		If[obj[[1]]==0 && obj[[2]]!=0,
			x1=xxmin;
			x2=xxmax;
			pd=Append[pd,plotdata[-obj[[3]]/obj[[2]],{x,x1,x2}]]
		];
		If[obj[[1]]!=0 && obj[[2]]!=0,
			xsect1= -(yymax*obj[[2]]+obj[[3]])/obj[[1]];
			xsect2= -(yymin*obj[[2]]+obj[[3]])/obj[[1]];
			x1=Max[xxmin,Min[xsect1,xsect2]];
			x2=Min[xxmax,Max[xsect1,xsect2]];
			pd=Append[pd,plotdata[-(obj[[1]]*xvar+obj[[3]])/obj[[2]],{xvar,x1,x2}]];
		];
	];
	If[args[[i]][[1]]=="Circle" && obj[[1]]!=0,(* if obj[[1]]==0 then it is a line.*) 
		option=args[[i]][[5]];
		r2=-obj[[4]]/obj[[1]]+(obj[[2]]^2+obj[[3]]^2)/(4*obj[[1]]^2);
		If[r2>0, pd=Append[pd,circledata[{-obj[[2]]/(2*obj[[1]]),-obj[[3]]/(2*obj[[1]])},Sqrt[r2]]]];
	];
	If[args[[i]][[1]]=="Mark",
		option=args[[i]][[5]];
		ptname="";ptposition=1;
		For[j=1,j<=Length[option],j++,
			Switch[option[[j]][[1]], 
				markname, ptname=option[[j]][[2]],
				markposition,ptposition=option[[j]][[2]]
			];
		];
		pd=Append[pd,anglemark[{obj[[1]][[1]]/obj[[1]][[3]],obj[[1]][[2]]/obj[[1]][[3]]},
		{obj[[2]][[1]]/obj[[2]][[3]],obj[[2]][[2]]/obj[[2]][[3]]},
		{obj[[3]][[1]]/obj[[3]][[3]],obj[[3]][[2]]/obj[[3]][[3]]}]];
		If[ptname!="", textdata=Append[textdata,Join[AngleMarkPosotion[{obj[[1]][[1]]/obj[[1]][[3]],obj[[1]][[2]]/obj[[1]][[3]]},
		{obj[[2]][[1]]/obj[[2]][[3]],obj[[2]][[2]]/obj[[2]][[3]]},
		{obj[[3]][[1]]/obj[[3]][[3]],obj[[3]][[2]]/obj[[3]][[3]]},ptposition],{"",ptname}]]];
	];
];
(* Print[windisp[pd]];*)
{pd,pointdata,textdata}
]

DrawConstruction[argseq__] := Module[{args,pd,pointdata,textdata}, 
args={argseq};
GetXYdomain[argseq];(*\:30aa\:30d7\:30b7\:30e7\:30f3\:3067\:306f\:305a\:305b\:308b\:3088\:3046\:306b\:3059\:308b\:ff0e*)
pd=MakePlotData[argseq];
pointdata=pd[[2]];
textdata=pd[[3]];
pd=pd[[1]];
setwindow[{XMIN,XMAX},{YMIN,YMAX}];
windisp[pd];
openfile[MyTeXFileName];
beginpicture["1cm"];
drwline[pd];
For[i=1,i<=Length[textdata],i++,
	letter[{textdata[[i]][[1]],textdata[[i]][[2]]},textdata[[i]][[3]],textdata[[i]][[4]]]
];
For[i=1,i<=Length[pointdata],i++,
	setpt[pointdata[[i]][[3]]];
	drwpt[{pointdata[[i]][[1]],pointdata[[i]][[2]]}]
];
endpicture[0];
closefile[];
]





