(* ::Package:: *)

(* ::Section:: *)
(*Mario's Mathematica Package*)


(* ::Subsubsection:: *)
(*Installation instructions*)


(* ::Text:: *)
(*Put the file in C:\Program Files\Wolfram Research\Mathematica\8.0\AddOns\Packages, and type*)


(* ::Text:: *)
(*Needs["Mario`"]*)


(* ::Text:: *)
(*to load the functions, or add Needs["Mario`"] to your C:\ProgramData\Mathematica\Kernel\init.m file to have the package automatically load on startup.*)


(* ::Subsection:: *)
(*Defined symbols*)


(* ::Subsubsection:: *)
(*Unit conversion:*)


(* ::Text:: *)
(*Convert[old,new] - converts numbers from one unit to another (of the same dimensionality); imported from Units`*)
(*ConvertFull[old,new] - My version of Convert[]; somewhat more successful with complicated derived SI units*)
(*ConvertAll[old,new] - converts numbers from one unit to any other, using \[HBar], c, G, Subscript[k, B], and Subscript[k, C] (the last is the Coulomb constant) to make the conversions from length to time to temperature, etc.*)
(*ConvertAllHL[old,new] - converts numbers from one unit to any other, using \[HBar], c, G, Subscript[k, B], and Subscript[\[CurlyEpsilon], 0] (the vacuum permittivity) to make the conversions from length to time to temperature, etc. HL stands for Heaviside-Lorentz, so named because this is the nondimensionalization convention of Heaviside-Lorentz CGS (the regular convention is Gaussian CGS).*)


(* ::Text:: *)
(*SI[expr] - converts numbers to SI units; modified from Units`*)
(*CGS[expr] - converts numbers to Gaussian CGS units (note: Ampere is not a fundamental unit in CGS; they are nondimensionalized by setting Subscript[k, C] to 1); modified from Units`*)
(*HeavisideLorentz[expr] - converts numbers to Heaviside-Lorentz CGS units, where Subscript[\[CurlyEpsilon], 0] nondimensionalizes Ampere*)
(*Planck[expr] - converts numbers to Planck units, where \[HBar], c, G, Subscript[k, B], and Subscript[k, C] are all equal to 1.*)
(*PlanckHL[expr] - converts numbers to Planck units, using the Heaviside-Lorentz convention: \[HBar], c, G, Subscript[k, B], and Subscript[\[CurlyEpsilon], 0] are all equal to 1.*)


(* ::Text:: *)
(*Const[expr] - converts actual constant variables (like \[HBar]) to their Mathematica forms (as in PlanckConstantReduced), and converts quoted unit names (as in "m") in the expression to their Mathematica form (Meter). Also works with SI prefixes, but the prefix must go inside the quotes (so "\[Micro]m" gives Micro Meter).*)
(*Clean[expr] - The "opposite" of Const[expr]: converts units to strings for presentation.*)
(*CConvert[old,new], CConvertFull[old,new], CConvertAll[old,new], CConvertAllHL[old,new] - versions of the corresponding functions above where the inputs have Const[] applied, and the result has Clean[] applied*)


(* ::Subsubsection:: *)
(*Complex plotting*)


(* ::Text:: *)
(*CPlot[expr,{z,xmin + I ymin,xmax + I ymax}] - Produce a colored complex plot of expr[z], with real bounds {xmin,xmax} and imaginary bounds {ymin,ymax}. Hue indicates the phase, with red meaning positive real, yellow-green meaning positive imaginary, cyan meaning negative real, and purple meaning negative imaginary, and with the luminance related to the absolute value of the number. Full color is \[LeftBracketingBar]z\[RightBracketingBar]=1, with black meaning 0 and white meaning infinite (in any direction).*)


(* ::Text:: *)
(*Options:*)


(* ::Text:: *)
(*HueSmoothing (default value 1): controls the smoothing of the hues; lower numbers mean more smoothing, range 0 to 2 (set to 0 to turn off)*)
(*GraySmoothing (default value 1/2): controls the smoothing of the luminance; lower numbers mean more smoothing, range 0 to 1 (set to 0 to turn off)*)
(*ZeroSize (default value 4/3): controls the "size" of zeros and poles (range 0 to \[Infinity]) - at low settings almost all numbers give full color, at high settings only values very close to 1 give color*)


(* ::Text:: *)
(*Method -> "Density" (default setting): use vertex coloring on a matrix of points; other options follow DensityPlot*)
(*Method -> "Overlay": use two DensityPlots on top of each other; other options follow Graphics*)
(*Method -> "Raster": use a Raster with colored rectangles (note: no smoothing in this mode); other options follow Graphics*)
(*PlotPoints: set number of sample points in the x and y directions (only use one number; separate x and y settings only supported in Method -> "Density")*)


(* ::Subsection:: *)
(*Code*)


BeginPackage["Mario`",{"Units`","PhysicalConstants`"}];
HueSmoothing::usage="HueSmoothing is an option for CPlot and related functions that chooses the amount of smoothing for the coloration: 2 for most smoothing, 0 for no smoothing.";
GraySmoothing::usage="GraySmoothing is an option for CPlot and related functions that chooses the amount of smoothing for the coloration: 1 for most smoothing, 0 for no smoothing.";
ZeroSize::usage="ZeroSize is an option for CPlot and related functions that chooses the relative strength of zeros and poles (black and white) as compared to midrange numbers (full color).";
CPlot::usage="\!\(CPlot[\*\nStyleBox[\"f\", \"TI\"], {\*\nStyleBox[\"z\", \"TI\"], \*SubscriptBox[\nStyleBox[\"z\", \"TI\"], \nStyleBox[\"min\", \"TI\"]], \*SubscriptBox[\nStyleBox[\"z\", \"TI\"], \nStyleBox[\"max\", \"TI\"]]}]\) generates a colored complex plot of \!\(\*\nStyleBox[\"f\", \"TI\"]\) as a function of the real and imaginary parts of \!\(\*\nStyleBox[\"z\", \"TI\"]\) from \!\(Re[\*SubscriptBox[\nStyleBox[\"z\", \"TI\"], \nStyleBox[\"min\", \"TI\"]]]\) to \!\(Re[\*SubscriptBox[\nStyleBox[\"z\", \"TI\"], \nStyleBox[\"max\", \"TI\"]]]\) along the \!\(\*StyleBox[\"x\", \"TI\"]\)-axis and from \!\(Im[\*SubscriptBox[\nStyleBox[\"z\", \"TI\"], \nStyleBox[\"min\", \"TI\"]]]\) to \!\(Im[\*SubscriptBox[\nStyleBox[\"z\", \"TI\"], \nStyleBox[\"max\", \"TI\"]]]\) along the \!\(\*StyleBox[\"y\", \"TI\"]\)-axis.";
Units`CGS::usage="CGS[expr] converts expr to Gaussian CGS units (centimeter/gram/second), where the unit of charge is chosen so that Coulomb's constant is 1.";
HeavisideLorentz::usage="HeavisideLorentz[expr] converts expr to Heaviside-Lorentz CGS units, where the unit of charge is chosen so that \!\(\*SubscriptBox[\(\[Epsilon]\), \(0\)]\)=1.";
Planck::usage="Planck[expr] converts expr to Planck units, where \[HBar] = c = G = \!\(\*SubscriptBox[\(k\), \(B\)]\) = \!\(\*SubscriptBox[\(k\), \(C\)]\) = 1.";
PlanckHL::usage="PlanckHL[expr] converts expr to Heaviside-Lorentz Planck units, where \[HBar] = c = G = \!\(\*SubscriptBox[\(k\), \(B\)]\) = \!\(\*SubscriptBox[\(\[Epsilon]\), \(0\)]\) = 1.";
ConvertFull::usage="ConvertFull[old, new] converts old to a form involving the combination of units new (and may have more success than Convert[old, new]).";
ConvertAll::usage="ConvertAll[old, new] converts old to a form involving the combination of units new, with no restriction on the dimensionality; extra factors of \[HBar], c, G, \!\(\*SubscriptBox[\(k\), \(B\)]\), and \!\(\*SubscriptBox[\(k\), \(C\)]\) are added as necessary to make the units come out right.";
ConvertAllHL::usage="ConvertAllHL[old, new] converts old to a form involving the combination of units new, with no restriction on the dimensionality; extra factors of \[HBar], c, G, \!\(\*SubscriptBox[\(k\), \(B\)]\), and \!\(\*SubscriptBox[\(\[Epsilon]\), \(0\)]\) are added as necessary to make the units come out right.";
Const::usage="Const[expr] converts all standard mathematical symbols like c, \[HBar], G, etc. to their proper values (in SI), and recognizes quoted units as their proper meaning (using the symbols from Units`).";
Clean::usage="Clean[expr] converts all units to their symbols, as embedded strings in the expression.";
(Evaluate[Symbol["C"~~#]]::usage="C"~~#~~"[old,new] is equivalent to "~~#~~"[old,new], but applies Const[] to its arguments and Clean[] to the result."
  )&/@{"Convert","ConvertFull","ConvertAll","ConvertAllHL"};
$ExtraSymbols::usage="$ExtraSymbols is a list of replacements of the form <symbol> \[Rule] <value> (for example, \[HBar] \[Rule] PlanckConstantReduced), which is added to the list of recognized symbols of Const[].";

Units`JulianYear::usage="JulianYear is a unit of time.";
Units`Franklin::usage="Franklin is the CGS unit of charge.";
Units`Stattesla::usage="Stattesla is a unit of magnetic flux density.";
Begin["`Private`"];


$prefixes={#[[1]],Evaluate[#[[1]]]::symbol}&/@Units`Private`$SIPrefixes[[1]]/."mu"->"\[Micro]";
{$replunit,$extractunit}=Dispatch/@{Rule@@#&/@#,Rule@@Reverse[#]&/@#}&@
  {{Abampere,"abA"},{Abcoulomb,"abC"},{Abfarad,"abF"},{Abhenry,"abH"},{Abohm,"ab\[CapitalOmega]"},{Abvolt,"abV"},{Ampere,"A"},{Angstrom,"\[Angstrom]"},{AstronomicalUnit,"AU"},
   {Atmosphere,"atm"},{AtomicMassUnit,"u"},{Bar,"bar"},{Barn,"b"},{Barye,"Ba"},{Becquerel,"Bq"},{Biot,"Bi"},{Bit,"bit"},{BritishThermalUnit,"Btu"},{BTU,"Btu"},
   {Calorie,"cal"},{Candela,"cd"},{Celsius,"\[Degree]C"},{Centimeter,"cm"},{Coulomb,"C"},{Curie,"Ci"},{Dalton,"Da"},{Day,"d"},{Dyne,"dyn"},{ElectronVolt,"eV"},{Erg,"erg"},
   {Fahrenheit,"\[Degree]F"},{Farad,"F"},{FluidOunce,"fl oz"},{Foot,"ft"},{Franklin,"Fr"},{Gal,"gal"},{Gallon,"gal"},{Gauss,"G"},{Gram,"g"},{GramWeight,"gf"},{Gravity,"g"},{GrayDose,"Gy"},
   {Hectare,"ha"},{Henry,"H"},{Hertz,"Hz"},{Horsepower,"hp"},{Hour,"h"},{Hundredweight,"cwt"},{ImperialPint,"ipt"},{Inch,"in"},{InchMercury,"inHg"},{Joule,"J"},
   {Kelvin,"K"},{Kilogram,"kg"},{KilogramForce,"kgf"},{LightYear,"ly"},{Liter,"l"},{Lumen,"lm"},{Lux,"lx"},{Maxwell,"Mx"},{Meter,"m"},{Mile,"mi"},
   {MillimeterMercury,"mmHg"},{Minute,"min"},{Mole,"mol"},{Month,"mo"},{Newton,"N"},{Oersted,"Oe"},{Ohm,"\[CapitalOmega]"},{Ounce,"oz"},{Parsec,"pc"},{Pascal,"Pa"},
   {Percent,"%"},{Phot,"ph"},{Pint,"pt"},{Poise,"P"},{Poise,"P"},{Pound,"lb"},{Poundal,"pdl"},{PoundForce,"lbf"},{PoundsPerSquareInch,"psi"},{Quart,"qt"},{Rad,"rad"},
   {Radian,"rad"},{Rankine,"R"},{Roentgen,"R"},{Rontgen,"R"},{Rutherford,"rd"},{Rydberg,"Ry"},{Second,"s"},{Siemens,"S"},{Slug,"slug"},{SolarMass,"M\[CircleDot]"},{Statampere,"statA"},
   {Statcoulomb,"statC"},{Statfarad,"statF"},{Stathenry,"statH"},{Statohm,"stat\[CapitalOmega]"},{Stattesla,"statT"},{Statvolt,"statV"},{Steradian,"sr"},{Stere,"st"},{Stilb,"sb"},{Stokes,"St"},
   {Stokes,"St"},{Tablespoon,"tbsp"},{Teaspoon,"tsp"},{Tesla,"T"},{Tonne,"t"},{UKGallon,"impgal"},{UKPint,"ipt"},{Volt,"V"},{Watt,"W"},{Weber,"Wb"},{Week,"wk"},
   {XUnit,"xu"},{Yard,"yd"},{Year,"yr"}};
$replpref=Dispatch@Quiet[(x_String #[[1]]->#[[2]]~~x&/@$prefixes)];
Begin["Global`"];
Mario`Private`$replconsts=Flatten[(u\[Function](#->u[[2]]&/@Flatten[{u[[1]]}]))/@Block[{$Context="Global`"},
  {g->AccelerationDueToGravity,{Subscript[N, A],NA}->AvogadroConstant,{Subscript[a, 0],a0}->BohrRadius,{k,Subscript[k, B],kB}->BoltzmannConstant,
   {Subscript[m, d],md}->DeuteronMass,{Subscript[m, E],mE}->EarthMass,{Subscript[R, E],RE}->EarthRadius,e->ElectronCharge,{Subscript[g, e],ge}->ElectronGFactor,
   Subscript[m, e]->ElectronMass,\[Alpha]->FineStructureConstant,G->GravitationalConstant,{Subscript[\[CapitalPhi], 0],\[CapitalPhi]0}->MagneticFluxQuantum,R->MolarGasConstant,
   {Subscript[V, m],Vm}->MolarVolume,{Subscript[m, \[Mu]],m\[Mu]}->MuonMass,{Subscript[m, n],mn}->NeutronMass,h->PlanckConstant,\[HBar]->PlanckConstantReduced,
   Subscript[m, P]->PlanckMass,{Subscript[m, p],mp}->ProtonMass,{Subscript[R, \[Infinity]],R\[Infinity]}->RydbergConstant,Ry->PlanckConstant RydbergConstant SpeedOfLight,
   {Subscript[R, sun],Rsun}->SolarRadius,c->SpeedOfLight,{Subscript[c, s],cs}->SpeedOfSound,\[Sigma]->StefanConstant,{Subscript[\[Mu], 0],\[Mu]0}->VacuumPermeability,
   {Subscript[\[Epsilon], 0],Subscript[\[CurlyEpsilon], 0],\[Epsilon]0,\[CurlyEpsilon]0}->VacuumPermittivity,{Subscript[\[Theta], W],\[Theta]W}->ArcSin[Sqrt[WeakMixingAngle]]}]];
End[];
$replconsts1=Join[{x:Subscript[_, _]->x},Cases[$replconsts,Except[Subscript[_,_]->_]]];
$replconsts2=Cases[$replconsts,x:(Subscript[_,_]->_)];
killPrefixes[s_String]:=Fold[If[MatchQ[#,_String]&&StringMatchQ[#,#2[[2]]~~__],#2[[1]]StringDrop[#,StringLength[#2[[2]]]],#]&,StringReplace[s,"\[Mu]"->"\[Micro]"],$prefixes]/.$extractunit
$ExtraSymbols={};
Const[expr_,protect_:{}]:=expr/.$extractunit/.x_String:>killPrefixes[x]/.
                    Join[#->#&/@protect,{x:Subscript[_, _]->x},Cases[$ExtraSymbols,Except[Subscript[_,_]->_]]]/.
                    Join[#->#&/@protect,Cases[$ExtraSymbols,Subscript[_,_]->_]]/.
                    Join[#->#&/@protect,$replconsts1]/.Join[#->#&/@protect,$replconsts2]
Clean[expr_]:=expr/.$replunit/.$replpref


Units`JulianYear::usage="JulianYear is a unit of time.";
$addenda={Mole->AvogadroConstant Mole,PrintersPoint->100/7227 Inch,ImperialPint->UKPint,UKPint->568261250Nano Liter,WineBottle->Gallon/5,Last->640UKGallon,
  ImperialGallon->UKGallon,UKGallon ->454609*^-5 Liter,Cup->231/16 Inch^3,Jigger -> 3/2 Shot,Pony -> 1/2 Jigger,Peck->215042/400 Inch^3,Omer->9/20Peck,
  SolarMass -> 1.98892*^33 Gram,AssayTon->(Milli Gram Ton)/TroyOunce,Grain->Pound/7000,Carat->1/5Gram,TroyOunce->480 Grain,Pennyweight->TroyOunce/20,AMU->AvogadroConstant^-1 Gram/Mole,
  Ounce->45359237/1600000 Gram,Poundal->Pound Foot/Second^2,TonForce->AccelerationDueToGravity LongTon,PoundForce->AccelerationDueToGravity Pound,
  KilogramForce->AccelerationDueToGravity Kilogram,Atmosphere->101325Pascal,InchMercury->Inch MillimeterMercury/(Milli Meter),PSI->PoundForce/Inch^2,
  MillimeterMercury->135951*^-8KilogramForce Centimeter^-2,Torr->Atmosphere/760,Rydberg->RydbergConstant PlanckConstant SpeedOfLight,Reyn->PoundForce Second/Inch^2,Rad->GrayDose/100,
  Rontgen->258*^-6 Coulomb/Kilogram,BohrMagneton->(ElectronCharge PlanckConstantReduced)/(2ElectronMass),NuclearMagneton->(ElectronCharge PlanckConstantReduced)/(2ProtonMass),JulianYear->36525/100Day,
  Parsec->648000/\[Pi] AU,Slug->PoundForce Second^2/Foot,Horsepower->550Foot PoundForce/Second,Statampere->SpeedOfLight^-1 (Meter Ampere)/(10Second),Statcoulomb->Statampere Second,
  Statfarad->Statcoulomb/Statvolt,Statohm->Statvolt/Statampere,Stathenry->Statohm Second,Stattesla->Statvolt Second/Meter^2,Statvolt->10^-6 Second/Meter SpeedOfLight Volt,Franklin->Statcoulomb};
$addenda=Thread[$addenda[[All,1]]->($addenda[[All,2]]//.Join[$addenda,Units`Private`$ToIntermediate[[1]]]//.Units`Private`$SIPrefixes//.Units`Private`$ToSI)];
Units`Private`$ToIntermediate=Dispatch@Evaluate@Join[FilterRules[Units`Private`$ToIntermediate[[1]],Except[$addenda]],$addenda];
Units`Private`$ToSI=Dispatch@Evaluate[Join[FilterRules[FilterRules[Units`Private`$ToSI[[1]],Except[$addenda]],Except[#]],#]&
  @{XUnit->1.00207699*^-13Meter,NauticalMile->1852 Meter,LightYear->31557600 SpeedOfLight Second,AstronomicalUnit -> 1.49597870700*^11 Meter,
    Gravity->AccelerationDueToGravity,Calorie->41868*^-4Joule,BTU->105505585262*^-8 Joule,ElectronVolt->ElectronCharge Volt}];


SI[expr_]:=expr//. Units`Private`$ToIntermediate//. Units`Private`$SIPrefixes//. Units`Private`$ToSI//.Units`Private`$ToFundamental
CGS[expr_]:=SI[expr]/.Ampere->10SpeedOfLight (Gram^(1/2) Centimeter^(3/2))/(Meter Second)//. Units`Private`$ToCGS
HeavisideLorentz[expr_]:=SI[expr]/.Ampere->10/Sqrt[4\[Pi]] SpeedOfLight (Gram^(1/2) Centimeter^(3/2))/(Meter Second)//. Units`Private`$ToCGS
Begin["Global`"];
Mario`Private`$Nondimensional=SI[Const[#]]&/@{Meter->Sqrt[Meter^2/(\[HBar] G/c^3)],Kilogram->Sqrt[Kilogram^2/(\[HBar] c/G)],Second->Sqrt[Second^2/(\[HBar] G/c^5)],Kelvin->Sqrt[Kelvin^2/(\[HBar] c^5/G/k^2)]};
End[];
$MeaninglessUnits={Bit->1,Radian->1,Steradian->1};
Planck[expr_]:=SI@CGS[expr/.Candela->Watt/Steradian]//.$MeaninglessUnits//.$Nondimensional
PlanckHL[expr_]:=SI@HeavisideLorentz[expr/.Candela->Watt/Steradian]//.$MeaninglessUnits//.$Nondimensional
ConvertFull[old_,new_]:=Units`Convert[SI[old],new]
ConvertAll[old_,new_]:=Planck[old]/Planck[new] new
ConvertAllHL[old_,new_]:=PlanckHL[old]/PlanckHL[new] new
(Evaluate[#][old_,new_]:=Clean[#2[Const@old,Const@new]])&[Symbol["C"~~#],Symbol[#]]&/@{"Convert","ConvertFull","ConvertAll","ConvertAllHL"};


kernel=Compile[{a,x},Piecewise[{{0,3x<=2(1-a)},{-1+(3x)/2,3x>=2+a}},(-2+2a+3x)^3/(54a^2)]];
hue=Compile[{n,x},Block[{m=4Mod[x,1],a=Min[n,1]},Cos[\[Pi] x]^2 (n-a)+(1+a-n)Piecewise[{{1-kernel[a,m],m<=1},{kernel[a,2-m],m<=2},{kernel[a,-2+m],m<=3}},1-kernel[a,4-m]]]];
huept[a_,x_]:={hue[a,x],hue[a,x-1/3],hue[a,x-2/3]}
huergb[a_,h_,white_,black_]:=1-black+huept[a,h] (white+black-1)
curve=Compile[{\[Alpha],a,x},2(-1+\[Alpha]) Piecewise[{{0,x<=-a},{x,x>=a}},a^-2 (1/8+(a-x)/(16 a))(a+x)^3]
  +Piecewise[{{1+2x-2x \[Alpha],x>=0},{(1+x)^2,\[Alpha]==0}},1-3x+2/\[Alpha]^2-4/\[Alpha]+(2x)/\[Alpha]+(-(2/\[Alpha]^2)+4/\[Alpha]) Sqrt[1+2x \[Alpha]-x \[Alpha]^2]]];
map[x_]=With[{f=(x\[Function]4/3 x-1/3 x^2)},\!\(\*
TagBox[GridBox[{
{"\[Piecewise]", GridBox[{
{
RowBox[{
RowBox[{"f", "[", "x", "]"}], "-", "1"}], 
RowBox[{"x", "<=", "1"}]},
{
RowBox[{"1", "-", 
RowBox[{"f", "[", 
RowBox[{"1", "/", "x"}], "]"}]}], "True"}
},
AllowedDimensions->{2, Automatic},
Editable->True,
GridBoxAlignment->{"Columns" -> {{Left}}, "ColumnsIndexed" -> {}, "Rows" -> {{Baseline}}, "RowsIndexed" -> {}},
GridBoxItemSize->{"Columns" -> {{Automatic}}, "ColumnsIndexed" -> {}, "Rows" -> {{1.}}, "RowsIndexed" -> {}},
GridBoxSpacings->{"Columns" -> {Offset[0.27999999999999997`], {Offset[0.84]}, Offset[0.27999999999999997`]}, "ColumnsIndexed" -> {}, "Rows" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}, "RowsIndexed" -> {}},
Selectable->True]}
},
GridBoxAlignment->{"Columns" -> {{Left}}, "ColumnsIndexed" -> {}, "Rows" -> {{Baseline}}, "RowsIndexed" -> {}},
GridBoxItemSize->{"Columns" -> {{Automatic}}, "ColumnsIndexed" -> {}, "Rows" -> {{1.}}, "RowsIndexed" -> {}},
GridBoxSpacings->{"Columns" -> {Offset[0.27999999999999997`], {Offset[0.35]}, Offset[0.27999999999999997`]}, "ColumnsIndexed" -> {}, "Rows" -> {Offset[0.2], {Offset[0.4]}, Offset[0.2]}, "RowsIndexed" -> {}}],
"Piecewise",
DeleteWithContents->True,
Editable->False,
SelectWithContents->True,
Selectable->False]\)];


Options[CPlotOverlay]=With[{o={HueSmoothing->1,GraySmoothing->1/2,ZeroSize->4/7,
  PlotPoints->40,Exclusions->None,PlotRange->All}},Sort@Join[o,FilterRules[Options[DensityPlot],Except[o]]]];
CPlotOverlay[f_,{z_,min_,max_},opts:OptionsPattern[CPlotOverlay]]:=
  Module[{x,y,q},q=f/. z->x+I y;
    Show[(DensityPlot[#[[1]],{x,Re[min],Re[max]},{y,Im[min],Im[max]},
        Exclusions->OptionValue[Exclusions],PlotRange->OptionValue[PlotRange],
        Evaluate[FilterRules[{opts},Options[Plot]]],
        ColorFunctionScaling->False,ColorFunction->#[[2]]]&)
      /@ {{Arg[N@q]/(2 \[Pi]),RGBColor@@huept[OptionValue[HueSmoothing],#]&},
         {map[Abs[q]],r\[Function]With[{\[Alpha]=OptionValue[ZeroSize],a=OptionValue[GraySmoothing]},
           With[{X=curve[\[Alpha],a,r],Y=curve[\[Alpha],a,-r]},
             GrayLevel[If[X+Y>=2,0,(-1+Y)/(-2+X+Y)],2-X-Y]]]}}]]
Options[CPlotDensity]=With[{o={HueSmoothing->1,GraySmoothing->1/2,ZeroSize->4/7,
  PlotPoints->40,Frame->True}},Sort@Join[o,FilterRules[Options[Graphics],Except[o]]]];
CPlotDensity[f_,{z_,min_,max_},opts:OptionsPattern[CPlotDensity]]:=
  Module[{x0=Re[min],x1=Re[max],y0=Im[min],y1=Im[max],approx,t,p=OptionValue[PlotPoints],
   color=z\[Function]With[{\[Alpha]=OptionValue[ZeroSize],a=OptionValue[GraySmoothing],r=map[Abs[N@z]]},
     RGBColor@@huergb[OptionValue[HueSmoothing],1/(2\[Pi]) Arg[N@z],curve[\[Alpha],a,r],curve[\[Alpha],a,-r]]]},
   t=Flatten[Table[{{x,y},color[f/.z->x+I y]},{x,x0,x1,(x1-x0)/(p-1)},{y,y0,y1,(y1-y0)/(p-1)}],1];
   Graphics[GraphicsComplex[t[[All,1]],Polygon[Flatten[Table[{0,1,p+1,p}+i+p j,{i,1,p-1},{j,0,p-2}],1]],
     VertexColors->t[[All,2]]],Frame->OptionValue[Frame],Evaluate[FilterRules[{opts},Options[Graphics]]]]]
Options[CPlotRaster]=With[{o={HueSmoothing->1,GraySmoothing->1/2,ZeroSize->4/7,
  PlotPoints->40,Frame->True}},Sort@Join[o,FilterRules[Options[Graphics],Except[o]]]];
CPlotRaster[f_,{z_,min_,max_},opts:OptionsPattern[CPlotRaster]]:=
  Module[{x0=Re[min],x1=Re[max],y0=Im[min],y1=Im[max],approx,t,p=OptionValue[PlotPoints],
    color=z\[Function]With[{\[Alpha]=OptionValue[ZeroSize],a=OptionValue[GraySmoothing],r=map[Abs[z]]},
      RGBColor@@huergb[OptionValue[HueSmoothing],1/(2\[Pi]) Arg[N@z],curve[\[Alpha],a,r],curve[\[Alpha],a,-r]]]},
    Graphics[Raster[ArrayPlot[Table[color[f/.z->x+I y],
        {y,y0+(y1-y0)/(2p),y1,(y1-y0)/p},{x,x0+(x1-x0)/(2p),x1,(x1-x0)/p}],
      DataReversed->True][[1,1]],{{x0,y0},{x1,y1}}],Frame->OptionValue[Frame],
      Evaluate[FilterRules[{opts},Options[Graphics]]]]]
CPlot[f_,{z_,min_,max_},opts___]:=
  Switch[Method/.{opts,Method->"Density"},
    "Overlay",CPlotOverlay,
    "Density",CPlotDensity,
    "Raster",CPlotRaster][f,{z,min,max},opts]


End[];
EndPackage[]
