(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



(* ::Title:: *)
(*concaveQuadrilateralDegreeOneInterpolants*)


(* ::Subtitle:: *)
(*Author: Gautam Dasgupta,*)
(*gd18@columbia.edu*)


(* ::Subsubtitle:: *)
(*July 5, 2007*)


(* ::Section:: *)
(*code*)


BeginPackage["concaveQuadrilateralDegreeOneInterpolants`","Global`"]


concaveQuadrilateralDegreeOneInterpolants::usage="concaveQuadrilateralDegreeOneInterpolants[nodes,{x,y}], second node to be the origin:
nodes: {{x1,y1},{0,0},{x3,y3},{x4,y4}}
generates degree one interpolant for all  nodes in x-y coordinates for a concave quadrilateral described by nodes."


s2::usage="s2[nodes,{x,y}], second node to be the origin:
nodes: {{x1,y1},{0,0},{x3,y3},{x4,y4}}
generates degree one interpolant only for the second node in x-y coordinates for a concave region described by nodes."


equationOfAnEdge::usage="
equationOfAnEdge[nodes,{x,y},{i,j}]  yields equation in x,y coordinates of the line joining nodes number i and j in the intercept form such that the value at the origin is 1."


alongAnEdge::usage="alongAnEdge[f,{x, y},nodes,{i,j},\[Tau]] yields a function of 0< \[Tau] < 1, of a function f in x,y coordinates between i and j nodes.

alongAnEdge[f,{x, y,z},nodes,{2,j},\[Tau]] also substitutes for z." 


equationOfAnEdge[nodes_,{x_,y_},{i_,j_}]:=Module[{s},s=Det[{{x,y,1},Append[nodes[[i]],1],Append[nodes[[j]],1]}];
(* check the sign *)
s / (s/.{x->0,y->0})//Expand]


alongAnEdge[f_,{x_, y_},nodes_,{i_,j_},\[Tau]_:\[Tau]]:=Module[{\[Tau]Rule},
\[Tau]Rule=Thread[{x,y}->nodes[[i]]+ \[Tau] (nodes[[j]]-nodes[[i]])
];
f/.\[Tau]Rule//Expand//Together//Chop]



alongAnEdge[f_,{x_, y_,z_},nodes_,{2,j_},\[Tau]_:\[Tau]]:=Module[{s},
s=alongAnEdge[f,{x, y},nodes,{2,j},\[Tau]];

(s/.z->\[Tau]*(Sqrt[(nodes[[2]]-nodes[[j]]).(nodes[[2]]-nodes[[j]])]))//Expand//Together//Chop
]



Begin["`Private`"]


s2[nodes_,{x_, y_}]:=Module[{z, num2,deno, \[Alpha], \[Beta], \[Gamma], x4,y4,eq1,s },

(* numerator as the product of side34 and side 41 *)
num2=Times@@(equationOfAnEdge[nodes,{x,y},#]& /@{{3,4},{4,1}});
deno=1+\[Alpha] x + \[Beta] y + \[Gamma]  z;

(* denominator should pass through the
downward projected point of node4 *)

{x4,y4}=nodes[[4]];
eq1=0==deno/.{x->x4, y->y4, z->-Sqrt[x4^2 + y4^2]};
deno=deno/.(Solve[eq1,\[Gamma]]//Flatten)[[1]];
s=(num2/deno);


(* along side21 and side23 from 2 the profile should be (1-\[Tau])
*)

along21=Cancel[alongAnEdge[s,{x, y,z},nodes,{2,1},\[Tau]]/(1-\[Tau])];

along23=Cancel[alongAnEdge[s,{x, y,z},nodes,{2,3},\[Tau]]/(1-\[Tau])];

(* unknown \[Alpha] and \[Beta] can now be solved *)

(s/.Flatten[
Solve[Thread[Cancel[Rationalize[{along21,along23}]]==1],{\[Alpha], \[Beta]}]
])/.z->Sqrt[x^2+y^2]

]



concaveQuadrilateralDegreeOneInterpolants[nodes_,{x,y}]:=
Module[{s,s1,s3,s4,eqHomogeneity,eqXLinearity,eqYLinearity,solS134 },
s=s2[nodes,{x,y}];
eqHomogeneity=s+s1+s3+s4==1;
eqXLinearity={s1, s, s3, s4}. First[Transpose[nodes]]==x;
eqYLinearity={s1, s, s3, s4}. Last[Transpose[nodes]]==y;
solS134=Solve[{eqHomogeneity, eqXLinearity, eqYLinearity},{s1,s3,s4}]//Together//Flatten;
{s1, s, s3, s4}={s1,s,s3, s4}/.solS134

]


End[]


EndPackage[]


(* ::Section:: *)
(*Example*)


(* ::Subsection:: *)
(*case-1*)


(* ::Input:: *)
(*nodes={{-1,0},{0,0},{-.5,-.5},{1,.5}}*)


(* ::Input:: *)
(*s=s2[nodes,{x,y}]*)


(* ::Subsection:: *)
(*case-2*)


(* ::Input:: *)
(*nodes={{-1,2},{0,0},{-.5,-.5},{1,.5}}*)


(* ::Input:: *)
(*Graphics[Polygon[nodes]]*)


(* ::Input:: *)
(*concaveQuadrilateralDegreeOneInterpolants[nodes,{x,y}]*)
