(*
:Titulo / Title: RLSCF3D
*)
(*
Representacion y localizacion de singularidades y ceros de funciones reales 
en 3D.
*)
(*
Representation and location of singularities and zeros of real functions in 
3D
*)
(*
:Autor / Author: Fernando Gomez Lanza
*)
(*
:email:  fgomez@pie.xtec.es
*)
(*
:URL:  http://www.xtec.es/~fgomez
*)
(*
Version: 1.0
*)
(*
:Fecha / Date: September 1999
*)
(*
:Requisitos / Requirements: Mathematica 3.0   ,   InequalitySolve.m 
*)
(*
:Nota: El analisis y representaci\[OAcute]n de funciones necesita Mathematica 
3.0. Este paquete proporciona buenos resultados para las funciones reales 3D 
mas comunes, como, polinomicas o fracion algebraica (incluyendo Hermite, 
Laguerre y cualquier otro tipo de polinomio), raices, trigonometricas e 
inversas, hiperbolicas e inversas, exponenciales y logaritmicas, utilizadas 
tambien dentro de Which[...] (tener cuidado al definir Which). No usar 
expresiones Which[...]/Which[...]. Las funciones que presentan iteraciones 
como Tan[Tan[x]], y combinaciones trigonometricas, 1/(Tan[x]+Sin[x]-1), dan 
el mismo resultado incorrecto que Plot3D, al aparecer ecuaciones 
trascendentes. En este caso debemos proceder manualmente. Podemos trabajar 
con otras funciones no definidas a\[NTilde]adiendo sus propiedades a 
ListaCeros y listaInfinitos. Ver el cuaderno de documentacion para 
procedimiento manual y otras explicaciones.
*)
(*
:Note: The analysis and representation of functions needs Mathematica 3.0. 
This package give us good results for common 3D real functions, such as, 
polynomial or algebric fraction (including Hermite, Laguerre, and all other 
polynomial type), roots, trigonometric and inverses, hiperbolic and inverses, 
exponential and logaritmic functions, using Which[..] also (take care with 
Which definitions). Not use Which[..]/Which[...] expressions. Functions with 
iterations, such as, Tan[Tan[x]], and trigonometric combinations, 1/(Tan[x]+
Sin[x]-1), give the same incorrect result with Plot3D, because appear 
trascendental equations. In this case we need to procede manually. We can 
work with other non defined functions adding the properties to ListaCeros and 
ListaInfinitos. See documentation notebook for manual procedure and other 
explanations.
*)
BeginPackage["RLSCF3D`"];

ListaCeros::usage=
  "ListaCeros={...,{f,v,p},...}  \n\n Lista de funciones analizadas para \
encontrar ceros. f=head-function, v=valor que cumple f[v]=0, p=periodo de la \
funcion f (0 si no es periodica). Ejemplo {...,{Tan,0,Pi},...}. \n\n List of \
analysed functions to find zeros. f=head-function, v=value that verifies \
f[v]=0, p=function's period (0 if not period). Example {...,{Tan,0,Pi},...} \
."; 
  
ListaInfinitos::usage=
  "ListaInfinitos={...,{f,v,p},...}  \n\n Lista de funciones analizadas para \
encontrar infinitos. f=head-function, v=valor que cumple f[v]=0, p=periodo de \
la funcion f (0 si no es periodica). Ejemplo {...,{Tan,Pi/2,Pi},...}. \n\n \
List of analysed functions to find infinities. f=head-function, v=value that \
verifies f[v]=Infinity, p=function's period (0 if not period). Example \
{...,{Tan,Pi/2,Pi},...} .";

Ceros::usage=
  "Ceros[f[x],{x,xmin,xmax}] \n\n Calcula los ceros de la funcion f[x]  entre \
xmin y xmax. \n\n Calculate zeros of function f[x] from xmin to xmax.";

Infinitos::usage=
  "Infinitos[f[x],{x,xmin,xmax}] \n\n Calcula los infinitos que aparecen en \
la funcion  f[x] entre xmin y xmax. \n\n Calculate infinities of function \
f[x] from xmin to xmax.";

Argum::usage=
  "Argum[f,expr] \n\n Encuentra los argumentos de la funcion f[ ] contenidas \
en la expresion expr. \n\n Finds the arguments of function f[ ] in expr. \n\n \
Argum[Tan,3x+5Tan[4x/8]].";

ArguWhich::usage=
  "ArguWhich[expr] \n\n Encuentra los argumentos de la funcion Which en la \
expresion expr. \n\n Finds the arguments of function Which in expr.";
	
LSFun3D::usage=
  "LSFun3D[f[x,y],{{x,xmin,xmax},{y,ymin,ymax}}] \n\n Proporciona la lista de \
posibles singularidades de la funcion f[x,y] . \n\n Provide the list of \
probable singularities of function f[x,y] . \n\n \
LSFun[Tan[x+y],{{x,-5,5},{y,0,10}}]";

SoloR::usage=
  "SoloR[f] \n\n Transforma la funcion f para que considere solo los valores \
reales, no considera las raices imaginarias de x^(1/3), por ejemplo. \n\n \
Transform the function f to consider real values only, example: complex roots \
of x^(1/3) are not considered.";

PlotFR3D::usage=
  "PlotFR3D[f[x,y],{{x,xmin,xmax},{y,ymin,ymax}}] \n\n Representa la version \
real de la funcion f[x,y]  evitando las discontinuidades. \n\n Plots the real \
version of function f[x,y]  avoiding discontinuities.";
 
dinmin::usage=
  "Longitud m\[IAcute]nima del intervalo a que se aplica PlotFR3D. ListaRec \
No genera rect\[AAcute]ngulos de longitud inferior a esta. Complementa a \
nitmax. \n\n Minimum longitude of the interval to that it is applied \
PlotFR3D. ListaRec doesn't generate   
  rectangles of inferior longitude to this. It supplements to nitmax.";

nitmax::usage=
  "M\[AAcute]ximo n\[UAcute]mero de divisiones del rect\[AAcute]ngulo donde \
se representa   
  la funci\[OAcute]n. Complementa a dinmin. \n\n Maximum number of divisions \
of the rectangle where it is represented   
  the function. It supplements to dinmin.";

vkmax::usage=
  "El intervalo {vkmin,vkmax} determinan los periodos calculados en presencia \
de funciones peri\[OAcute]dicas. Ambos deben ser n\[UAcute]meros enteros. \n\n\
 The interval {vkmin,vkmax} determine the periods calculated in presence of   \

periodic functions. Both should be integer numbers."
	
vkmin::usage=
  "El intervalo {vkmin,vkmax} determinan los periodos calculados en presencia \
de funciones peri\[OAcute]dicas. Ambos deben ser n\[UAcute]meros enteros. \n\n\
 The interval {vkmin,vkmax} determine the periods calculated in presence of   \

periodic functions. Both should be integer numbers."

Sensibilidad::usage=
  "Determina el n\[UAcute]mero de puntos que verifican la presencia de la \
singularidad en una zona cuando estas aparecen como funciones no polinomicas. \
\n\n It determines the number of points where the program verifies the \
presence of the singularity in an area when these appear in form of non \
polynomial functions."

dinmin=1/10;
nitmax=10;
vkmax=10;vkmin=-10;
Sensibilidad=100;
		
ListaCeros={{Null,0,0},{Log,1,0},{Tan,0,Pi},{Cot,Pi/2,Pi},{Sin,0,Pi},{Cos,
      Pi/2,Pi},{Sinh,0,0},{Tanh,0,0}};

ListaInfinitos={{Log,0,0},{Tan,Pi/2,Pi},{Cot,0,Pi},{Csc,0,Pi},{Sec,Pi/2,Pi},{
      Csch,0,0},{Coth,0,0}};

Unprotect[Log];
Log[Global`x_/Global`y_]=Log[Global`x]-Log[Global`y];
Log[Global`b_,Global`x_/Global`y_]=
  Log[Global`b,Global`x]-Log[Global`b,Global`y];
Protect[Log];
 
Begin["`privado`"]; 

Argum[f_,expr_]:=
  If[f===Null,{expr},
    Union[Map[Apply[Part,Flatten[{expr,#,1}]]&,Position[expr,f[_]]]]];

ArguWhich[f_]:=
  Flatten[Map[Partition[Apply[List,Apply[Part,#]],2]&,
      Map[Drop[Prepend[#,f],-1]&,Position[f,Which]]],1];

DescomponerWhich[f_/;FreeQ[f,Which]==False]:=
  ReplaceAll[f,Which[__,___]->EliRepes[Map[Last[#]&,ArguWhich[f]]]];

ArguPower[f_]:=
  Flatten[Map[Apply[List,#]&,
      DeleteCases[
        Map[Apply[Part,#]&,Map[Drop[Prepend[#,f],-1]&,Position[f,Power]]],
        Power[_,_Integer]]],1];

Needs["Algebra`InequalitySolve`"];

EliRepes[l_List]:=Union[l,l];

LSGen3D[ff_,l_List:{{Global`x,-10,10}},g_List:{Null,0,0}]:=
  Module[{p,s,lx},Clear[k];lx=Map[First[#]&,l];p=Argum[g[[1]],ff];
    s=Flatten[Map[Solve[#==g[[2]]+k*g[[3]],lx]&,p]];
		s=If[FreeQ[s,k],s,Flatten[Table[s,{k,vkmin,vkmax}],1]];s=Map[(lx/.#)&,s];
    s=Chop[N[s],10^(-$MachinePrecision+1)];s];

Ceros[f_,l_List]:=EliRepes[Flatten[Map[LSGen3D[f,l,#]&,ListaCeros],1]];

Infinitos[f_,l_List]:=
  EliRepes[Flatten[Map[LSGen3D[f,l,#]&,ListaInfinitos],1]];

LSFun3D[f_,l_List:{{Global`x,-10,10},{Global`y,-10,10}}]:=
  Module[{lsum,ls},lsum=Apart[f];
    lsum=Switch[lsum,_Plus,Apply[List,lsum],_,{lsum}];
			ls=If[FreeQ[f,Which],{},
        Map[Map[First,l]/.#&,
          EliRepes[
            Map[InequalitySolve[#,Map[First,l]]&,
                Flatten[
                  Map[First,ArguWhich[f]]/.{
                      Less[a_,b_,c__]->{Less[a,b],Less[b,c]},
                      Greater[a_,b_,c__]->{Greater[a,b],Greater[b,c]},
                      Equal[a_,b_,c__]->{Equal[a,b],Equal[b,c]},
                      And->List}]]/.{Less->Rule,Greater->Rule,
                Equal->Rule}]]];
		lsum=Flatten[Map[If[FreeQ[#,Which],#,DescomponerWhich[#]]&,lsum]];
    lsum=Union[lsum,Flatten[Map[ArguPower[#]&,lsum],1]];
		EliRepes[
      Flatten[{{ls},Map[Ceros[#,l]&,EliRepes[Map[Denominator[#]&,lsum]]],
          Map[Infinitos[#,l]&,EliRepes[Map[Numerator[#]&,lsum]]]},2]]];

NRcs[ss_:{Null,Null},lxy_]:=
  Module[{i,v,r=False},v=If[FreeQ[ss,lxy[[2,1]]],lxy[[1]],lxy[[2]]];
    For[i=0,(i<Sensibilidad)&&(r==False),i++,
      r=If[((Im[ss]=={0,0})&&(lxy[[1,2]]<ss[[1]]<lxy[[1,3]])&&(
                  lxy[[2,2]]<ss[[2]]<lxy[[2,3]]))/.v[[1]]->
              N[i*(v[[3]]-v[[2]])/Sensibilidad+v[[2]]],True,False,False]];r];

Rcs[ss_:{Null,Null},lxy_]:=
  If[PolynomialQ[ss[[1]],{lxy[[1,1]],lxy[[2,1]]}]&&
      PolynomialQ[ss[[2]],{lxy[[1,1]],lxy[[2,1]]}],
    InequalitySolve[(lxy[[1,2]]<ss[[1]]<lxy[[1,3]])&&(
          lxy[[2,2]]<ss[[2]]<lxy[[2,3]]),
      If[FreeQ[ss,lxy[[2,1]]],lxy[[1,1]],lxy[[2,1]]]],NRcs[ss,lxy]];

DivRect[r_:{{-1,1},{-1,1}},dmin_,n_Integer:2]:=
  Module[{dx,dy,nx,ny},dx=Apply[Subtract,{r[[1,2]],r[[1,1]]}]/n;
    dy=Apply[Subtract,{r[[2,2]],r[[2,1]]}]/n;If[dx<dmin,nx=1;dx=dx*n;,nx=n];
    If[dy<dmin,ny=1;dy=dy*n;,ny=n];
    Flatten[Table[{{r[[1,1]]+(i-1)*dx,r[[1,1]]+i*dx},{r[[2,1]]+(j-1)*dy,
            r[[2,1]]+j*dy}},{i,1,nx},{j,1,ny}],1]];

ListaRec[ls_,lxy_:{{x,-10,10},{y,-10,10}},dmin_:1/10,nmax_:10]:=
  Module[{q,nls,n=0,lr={},lrq={},lrq0={}},lrq={Map[Rest,lxy]};nls=Length[ls];
		While[(Length[lrq]>0)&&(n<nmax),lrq0=lrq;lrq={};n++;
			For[i=1,i<=Length[lrq0],
				q=Count[
            Map[Rcs[#,{Prepend[lrq0[[i,1]],lxy[[1,1]]],
                    Prepend[lrq0[[i,2]],lxy[[2,1]]]}]&,ls],False];
        Switch[q,nls,lr=Union[lr,{lrq0[[i]]}],_,
          lrq=Union[lrq,DivRect[lrq0[[i]],dmin,2]]];i++]
		];Map[{Join[{lxy[[1,1]]},#[[1]]],Join[{lxy[[2,1]]},#[[2]]]}&,lr]
	];

PlotPuntos[lxy0_,lxy_,p0_:15]:=
  IntegerPart[
    N[{3,3}+(p0-3)*{(lxy[[1,3]]-lxy[[1,2]])/(lxy0[[1,3]]-lxy0[[1,2]]),(
                lxy[[2,3]]-lxy[[2,2]])/(lxy0[[2,3]]-lxy0[[2,2]])}]];

SoloR[f_]:=
  f//.a_^n_:>
      ((Sign[a])^Numerator[Rationalize[n]])*Abs[a]^
            n/;((Head[a]=!=Abs)&&(Mod[Denominator[Rationalize[n]],2]=!=0)&&(
              IntegerQ[n]==False));

PlotFR3D[f_,l_List:{{Global`x,-10,10},{Global`y,-10,10}},opciones___]:=
  Module[{ff,lsin,lrec,dibn},SetOptions[Plot3D,opciones];ff=SoloR[f];
    lsin=LSFun3D[ff,l];lrec=ListaRec[lsin,l,dinmin,nitmax];
    dibn=Map[Plot3D[ff,Evaluate[First[#]],Evaluate[Last[#]],
            PlotPoints->
              Evaluate[PlotPuntos[l,#,Options[Plot3D,PlotPoints][[1,2]]]],
            DisplayFunction->Identity]&,lrec];
    Show[dibn,DisplayFunction->$DisplayFunction]];

End[];

EndPackage[];
