Index: trunk/examples/extended/electromagnetic/TestEm7/include/DetectorConstruction.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/DetectorConstruction.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/DetectorConstruction.hh	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: DetectorConstruction.hh,v 1.2 2006/06/29 16:57:29 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/DetectorMessenger.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/DetectorMessenger.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/DetectorMessenger.hh	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: DetectorMessenger.hh,v 1.3 2006/06/29 16:57:31 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/EventAction.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/EventAction.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/EventAction.hh	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: EventAction.hh,v 1.2 2006/06/29 16:57:34 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/EventActionMessenger.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/EventActionMessenger.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/EventActionMessenger.hh	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: EventActionMessenger.hh,v 1.3 2006/06/29 16:57:36 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/G4ScreenedNuclearRecoil.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/G4ScreenedNuclearRecoil.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/G4ScreenedNuclearRecoil.hh	(revision 1230)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4ScreenedNuclearRecoil.hh,v 1.3 2007/12/07 17:51:10 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// G4ScreenedNuclearRecoil.hh,v 1.24 2008/05/01 19:58:59 marcus Exp
+// GEANT4 tag 
 //
 //
@@ -45,4 +45,5 @@
 //
 // First version, April 2004, Marcus H. Mendenhall, Vanderbilt University
+// May 1, 2008 -- Added code to allow process to have zero cross section above max energy, to coordinate with G4MSC.  -- mhm
 //
 // Class Description - End
@@ -54,4 +55,5 @@
 #include "globals.hh"
 #include "G4VDiscreteProcess.hh"
+#include "G4ParticleChange.hh"
 #include "c2_function.hh"
 
@@ -59,9 +61,13 @@
 #include <vector>
 
-class G4VParticleChange;
+class G4VNIELPartition;
+
+typedef c2_const_ptr<G4double> G4_c2_const_ptr;
+typedef c2_ptr<G4double> G4_c2_ptr;
+typedef c2_function<G4double> G4_c2_function;
 
 typedef struct G4ScreeningTables { 
 	G4double z1, z2, m1, m2, au, emin;
-	c2_function<G4double> *EMphiData; 
+	G4_c2_const_ptr EMphiData; 
 } G4ScreeningTables;
 
@@ -70,12 +76,11 @@
 {
 public:
-        G4ScreenedCoulombCrossSectionInfo() { }
-        ~G4ScreenedCoulombCrossSectionInfo() { }
-
-        const char *CVSHeaderVers() { return
-                "";
-        }
-
-        const char *CVSFileVers() { return ""; }
+	G4ScreenedCoulombCrossSectionInfo() { }
+	~G4ScreenedCoulombCrossSectionInfo() { }
+	
+	static const char* CVSHeaderVers() { return 
+		"G4ScreenedNuclearRecoil.hh,v 1.24 2008/05/01 19:58:59 marcus Exp GEANT4 tag ";
+	}
+	static const char* CVSFileVers();
 };
 
@@ -114,6 +119,6 @@
 	
 	// get the mean-free-path table for the indexed material 
-	c2_function<G4double> * operator [] (G4int materialIndex) { 
-			return MFPTables.find(materialIndex)!=MFPTables.end() ? MFPTables[materialIndex] : (c2_function<G4double> *)0;
+	const G4_c2_function * operator [] (G4int materialIndex) { 
+			return MFPTables.find(materialIndex)!=MFPTables.end() ? &(MFPTables[materialIndex].get()) : (G4_c2_function *)0;
 		}
 	
@@ -122,6 +127,6 @@
 	ParticleCache targetMap;
 	G4int verbosity;
-	std::map<G4int, c2_function<G4double> *> sigmaMap; // total cross section for each element
-	std::map<G4int, c2_function<G4double> *> MFPTables; // MFP for each material
+	std::map<G4int, G4_c2_const_ptr > sigmaMap; // total cross section for each element
+	std::map<G4int, G4_c2_const_ptr > MFPTables; // MFP for each material
 	
 private:
@@ -134,5 +139,7 @@
 	G4ScreenedCoulombCrossSection *crossSection;
 	G4double a1, a2, sinTheta, cosTheta, sinZeta, cosZeta, eRecoil;
-	G4ParticleDefinition *recoilIon;	} G4CoulombKinematicsInfo;
+	G4ParticleDefinition *recoilIon;	
+	const G4Material *targetMaterial;
+} G4CoulombKinematicsInfo;
 
 class G4ScreenedCollisionStage {
@@ -146,5 +153,5 @@
 
 public:
-	G4ScreenedCoulombClassicalKinematics() { }
+	G4ScreenedCoulombClassicalKinematics();
 	virtual void DoCollisionStep(class G4ScreenedNuclearRecoil *master,
 		const class G4Track& aTrack, const class G4Step& aStep);
@@ -153,5 +160,11 @@
 		const G4ScreeningTables *screen, 
 		G4double eps, G4double beta);
-	virtual ~G4ScreenedCoulombClassicalKinematics() {}
+	virtual ~G4ScreenedCoulombClassicalKinematics() { }
+protected:
+	// the c2_functions we need to do the work.
+	c2_const_plugin_function_p<G4double> &phifunc;
+	c2_linear_p<G4double> &xovereps;
+	G4_c2_ptr diff;
+	
 };
 
@@ -165,4 +178,9 @@
 };
 
+/**
+	\brief A process which handles screened Coulomb collisions between nuclei
+ 
+*/
+
 class G4ScreenedNuclearRecoil : public G4ScreenedCoulombCrossSectionInfo, public G4VDiscreteProcess
 {
@@ -171,54 +189,120 @@
 	friend class G4ScreenedCollisionStage;
 
+	/// \brief Construct the process and set some physics parameters for it.
+	/// \param processName the name to assign the process
+	/// \param ScreeningKey the name of a screening function to use.  
+	/// The default functions are "zbl" (recommended for soft scattering),
+	/// "lj" (recommended for backscattering) and "mol" (Moliere potential)
+	/// \param GenerateRecoils if frue, ions struck by primary are converted into new moving particles.
+	/// If false, energy is deposited, but no new moving ions are created.
+	/// \param RecoilCutoff energy below which no new moving particles will be created, 
+	/// even if \a GenerateRecoils is true.
+	/// Also, a moving primary particle will be stopped if its energy falls below this limit.
+	/// \param PhysicsCutoff the energy transfer to which screening tables are calucalted.  
+	/// There is no really
+	/// compelling reason to change it from the 10.0 eV default.  However, see the paper on running this
+	/// in thin targets for further discussion, and its interaction with SetMFPScaling()
 	G4ScreenedNuclearRecoil(const G4String& processName = "ScreenedElastic",
 			const G4String &ScreeningKey="zbl", G4bool GenerateRecoils=1, 
 			G4double RecoilCutoff=100.0*eV, G4double PhysicsCutoff=10.0*eV);
-	
+	/// \brief destructor
 	virtual ~G4ScreenedNuclearRecoil();
-	
+	/// \brief used internally by Geant4 machinery
 	virtual G4double GetMeanFreePath(const G4Track&, G4double, G4ForceCondition* );
-	
+	/// \brief used internally by Geant4 machinery
 	virtual G4VParticleChange* PostStepDoIt(const G4Track& aTrack, const G4Step& aStep);
-	
+	/// \brief test if a prticle of type \a aParticleType can use this process
+	/// \param aParticleType the particle to test
 	virtual G4bool IsApplicable(const G4ParticleDefinition& aParticleType);
-	
+	/// \brief Build physics tables in advance.  Not Implemented.
+	/// \param aParticleType the type of particle to build tables for 
 	virtual void BuildPhysicsTable(const G4ParticleDefinition&) { }
-	
+	/// \brief Export physics tables for persistency.  Not Implemented.
+	/// \param aParticleType the type of particle to build tables for 
 	virtual void DumpPhysicsTable(const G4ParticleDefinition& aParticleType);
-	
+	/// \brief deterine if the moving particle is within  the strong force range of the selected nucleus
+	/// \param A the nucleon number of the beam
+	/// \param A1 the nucleon number of the target
+	/// \param apsis the distance of closest approach
 	virtual G4bool CheckNuclearCollision(G4double A, G4double A1, G4double apsis); // return true if hard collision
 
 	virtual G4ScreenedCoulombCrossSection *GetNewCrossSectionHandler(void);
 	
-	G4double GetNIEL() const { return NIEL; } // Get non-ionizing energy loss for last step
-	
+	/// \brief Get non-ionizing energy loss for last step
+	G4double GetNIEL() const { return NIEL; } 
+	
+	/// \brief clear precomputed screening tables
 	void ResetTables(); // clear all data tables to allow changing energy cutoff, materials, etc.
-	
+
+	/// \brief set the upper energy beyond which this process has no cross section
+	///
+	/// This funciton is used to coordinate this process with G4MSC.  Typically, G4MSC should 
+	///  not be allowed to operate in a range which overlaps that of this process.  The criterion which is most reasonable
+	/// is that the transition should be somewhere in the modestly relativistic regime (500 MeV/u for example).
+	/// \param energy energy per nucleon for the cutoff
+	void SetMaxEnergyForScattering(G4double energy) { processMaxEnergy=energy; }
+	/// \brief find out what screening funciton we are using
 	std::string GetScreeningKey() const { return screeningKey; }
+	/// \brief enable or disable all energy deposition by this process
+	/// \param flag if true, enable deposition of energy (the default).  If false, disable deposition.
 	void AllowEnergyDeposition(G4bool flag) { registerDepositedEnergy=flag; }
+	/// \brief get flag indicating whether deposition is enabled
 	G4bool GetAllowEnergyDeposition() const { return registerDepositedEnergy; }
+	/// \brief enable or disable the generation of recoils.  
+	/// If recoils are disabled, the energy they would have received is just deposited.
+	/// \param flag if true, create recoil ions in cases in which the energy is above the recoilCutoff.  
+	/// If false, just deposit the energy.
 	void EnableRecoils(G4bool flag) { generateRecoils=flag; }
+	/// \brief find out if generation of recoils is enabled.
 	G4bool GetEnableRecoils() const { return generateRecoils; }
+	/// \brief set the mean free path scaling as specified
+	/// \param scale the factor by which the default MFP will be scaled.  
+	/// Set to less than 1 for very thin films, typically, to sample multiple scattering,
+	/// or to greater than 1 for quick simulaitons with a very long flight path.
 	void SetMFPScaling(G4double scale) { MFPScale=scale; }
+	/// \brief get the MFPScaling parameter
 	G4double GetMFPScaling() const { return MFPScale; }
+	/// \brief enable or disable whether this process will skip collisions 
+	/// which are close enough they need hadronic phsyics. Default is true (skip close collisions).
+	/// Disabling this results in excess nuclear stopping power.
+	/// \param flag true results in hard collisions being skipped.  false allows hard collisions.
 	void AvoidNuclearReactions(G4bool flag) { avoidReactions=flag; }
+	/// \brief get the flag indicating whether hadronic collisions are ignored.
 	G4bool GetAvoidNuclearReactions() const { return avoidReactions; }
+	/// \brief set the minimum energy (per nucleon) at which recoils can be generated, 
+	/// and the energy (per nucleon) below which all ions are stopped. 
+	/// \param energy energy per nucleon 
 	void SetRecoilCutoff(G4double energy) { recoilCutoff=energy; }
+	/// \brief get the recoil cutoff
 	G4double GetRecoilCutoff() const { return recoilCutoff; }
+	/// \brief set the energy to which screening tables are computed.  Typically, this is 10 eV or so, and not often changed.
+	/// \param energy the cutoff energy 
 	void SetPhysicsCutoff(G4double energy) { physicsCutoff=energy; ResetTables(); }
+	/// \brief get the physics cutoff energy.
 	G4double GetPhysicsCutoff() const { return physicsCutoff; }
-	class G4ParticleChange &GetParticleChange() { return aParticleChange; }
-	void AddToNIEL(G4double energy) { NIEL+=energy; }
+	/// \brief set the pointer to a class for paritioning energy into NIEL
+	/// \brief part the pointer to the class.
+	void SetNIELPartitionFunction(const G4VNIELPartition *part);
+	/// \brief set the cross section boost to provide faster computation of backscattering
+	/// \param fraction the fraction of particles to have their cross section boosted.
+	/// \param HardeningFactor the factor by which to boost the scattering cross section.
 	void SetCrossSectionHardening(G4double fraction, G4double HardeningFactor) {
 		hardeningFraction=fraction;
 		hardeningFactor=HardeningFactor;
 	}
+	/// \brief get the fraction of particles which will have boosted scattering
 	G4double GetHardeningFraction() const { return hardeningFraction; }
+	/// \brief get the boost factor in use.
 	G4double GetHardeningFactor() const { return hardeningFactor; }
+	/// \brief the the interaciton length used in the last scattering.
 	G4double GetCurrentInteractionLength() const { return currentInteractionLength; }
+	/// \brief set a function to compute screening tables, if the user needs non-standard behavior.
+	/// \param cs a class which constructs the screening tables.
 	void SetExternalCrossSectionHandler(G4ScreenedCoulombCrossSection *cs) { 
 		externalCrossSectionConstructor=cs; 
 	}
+	/// \brief get the verbosity.
 	G4int GetVerboseLevel() const { return verboseLevel; }
+
 	std::map<G4int, G4ScreenedCoulombCrossSection*> &GetCrossSectionHandlers() 
 		{ return crossSectionHandlers; }
@@ -228,12 +312,22 @@
 	void SetValidCollision(G4bool flag) { validCollision=flag; }
 	G4bool GetValidCollision() const { return validCollision; }
+
+	/// \brief get the pointer to our ParticleChange object.  for internal use, primarily.
+	class G4ParticleChange &GetParticleChange() { return static_cast<G4ParticleChange &>(*pParticleChange); }
+	/// \brief take the given energy, and use the material information to partition it into NIEL and ionizing energy.
+	void DepositEnergy(G4int z1, G4double a1, const G4Material *material, G4double energy);
 	
 protected:
-	G4double highEnergyLimit, lowEnergyLimit;
+	/// \brief the energy per nucleon above which the MFP is constant 
+	G4double highEnergyLimit;
+	/// \brief the energy per nucleon below which the MFP is zero
+	G4double lowEnergyLimit;
+	/// \brief the energy per nucleon beyond which the cross section is zero, to cross over to G4MSC
+	G4double processMaxEnergy;
 	G4String screeningKey;
 	G4bool generateRecoils, avoidReactions;
 	G4double recoilCutoff, physicsCutoff;
 	G4bool registerDepositedEnergy;
-	G4double NIEL;
+	G4double IonizingLoss, NIEL;
 	G4double MFPScale;
 	G4double hardeningFraction, hardeningFactor;
@@ -243,9 +337,8 @@
 	
 	std::map<G4int, G4ScreenedCoulombCrossSection*> crossSectionHandlers;		
-	std::map<G4int, c2_function<G4double>*> meanFreePathTables;
 	
 	G4bool validCollision;
 	G4CoulombKinematicsInfo kinematics;
-
+	const G4VNIELPartition *NIELPartitionFunction;
 };
 
@@ -268,5 +361,5 @@
 	std::vector<G4String> GetScreeningKeys() const;
 	
-	typedef c2_function<G4double> &(*ScreeningFunc)(G4int z1, G4int z2, size_t nPoints, G4double rMax, G4double *au);
+	typedef G4_c2_function &(*ScreeningFunc)(G4int z1, G4int z2, size_t nPoints, G4double rMax, G4double *au);
 	
 	void AddScreeningFunction(G4String name, ScreeningFunc fn) {
@@ -279,5 +372,8 @@
 };
 
-
+G4_c2_function &ZBLScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval);
+G4_c2_function &MoliereScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval);
+G4_c2_function &LJScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval);
+G4_c2_function &LJZBLScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double *auval);
 
 #endif
Index: trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmLivermore.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmLivermore.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmLivermore.hh	(revision 1230)
@@ -26,5 +26,5 @@
 //
 // $Id: PhysListEmLivermore.hh,v 1.1 2006/11/22 18:56:20 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmPenelope.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmPenelope.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmPenelope.hh	(revision 1230)
@@ -26,5 +26,5 @@
 //
 // $Id: PhysListEmPenelope.hh,v 1.1 2006/11/22 18:56:21 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmStandard.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmStandard.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmStandard.hh	(revision 1230)
@@ -26,5 +26,5 @@
 //
 // $Id: PhysListEmStandard.hh,v 1.3 2006/06/29 16:57:44 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmStandardNR.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmStandardNR.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmStandardNR.hh	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: PhysListEmStandardNR.hh,v 1.1 2008/01/14 12:11:38 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmStandardSS.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmStandardSS.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmStandardSS.hh	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: PhysListEmStandardSS.hh,v 1.1 2006/10/24 11:37:56 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/PhysicsList.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/PhysicsList.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/PhysicsList.hh	(revision 1230)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: PhysicsList.hh,v 1.6 2006/11/22 18:56:21 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// $Id: PhysicsList.hh,v 1.8 2008/11/20 20:34:50 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -39,4 +39,5 @@
 
 #include "G4VModularPhysicsList.hh"
+#include "G4EmConfigurator.hh"
 #include "globals.hh"
 
@@ -50,4 +51,5 @@
 {
 public:
+
   PhysicsList();
   virtual ~PhysicsList();
@@ -67,4 +69,7 @@
 
 private:
+
+  G4EmConfigurator em_config; 
+
   G4double cutForGamma;
   G4double cutForElectron;
@@ -77,4 +82,5 @@
   G4String                             emName;
   G4VPhysicsConstructor*               emPhysicsList;
+  G4VPhysicsConstructor*               decPhysicsList;
   std::vector<G4VPhysicsConstructor*>  hadronPhys;
     
Index: trunk/examples/extended/electromagnetic/TestEm7/include/PhysicsListMessenger.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/PhysicsListMessenger.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/PhysicsListMessenger.hh	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: PhysicsListMessenger.hh,v 1.3 2006/06/29 16:57:52 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/PrimaryGeneratorAction.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/PrimaryGeneratorAction.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/PrimaryGeneratorAction.hh	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: PrimaryGeneratorAction.hh,v 1.3 2006/06/29 16:57:54 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/PrimaryGeneratorMessenger.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/PrimaryGeneratorMessenger.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/PrimaryGeneratorMessenger.hh	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: PrimaryGeneratorMessenger.hh,v 1.3 2006/06/29 16:57:56 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/RunAction.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/RunAction.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/RunAction.hh	(revision 1230)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: RunAction.hh,v 1.12 2008/01/14 12:11:39 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: RunAction.hh,v 1.14 2008/08/22 18:30:27 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -40,10 +40,5 @@
 class PrimaryGeneratorAction;
 class G4Run;
-
-namespace AIDA {
- class IAnalysisFactory;
- class ITree;
- class IHistogram1D;
-} 
+class Histo;
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -52,5 +47,6 @@
 {
 public:
-  RunAction(DetectorConstruction*, PhysicsList*,PrimaryGeneratorAction*);
+  RunAction(DetectorConstruction*, PhysicsList*,
+	    PrimaryGeneratorAction*);
   virtual ~RunAction();
 
@@ -58,20 +54,19 @@
   void   EndOfRunAction(const G4Run*);
     
-  void FillTallyEdep(G4int n, G4double e) {tallyEdep[n] += e;};
+  void FillTallyEdep(G4int n, G4double e)  {tallyEdep[n] += e;};
   void FillEdep(G4double de, G4double eni) {edeptot += de; eniel += eni;};
        
   G4double GetBinLength() {return binLength;};
   G4double GetLength()    {return length;};
-  G4double GetOffsetX()   {return offsetX;} 
-  void     FillHisto(G4int id, G4double x, G4double weight = 1.0);
+  G4double GetOffsetX()   {return offsetX;}
+ 
+  void FillHisto(G4int id, G4double x, G4double weight = 1.0);
     
-  void AddProjRange (G4double x) {projRange += x; projRange2 += x*x;};
+  void AddProjRange (G4double x) 
+  {projRange += x; projRange2 += x*x; nRange++;};
   void AddPrimaryStep() {nPrimarySteps++;};
                    
 private:  
-  void bookHisto();
-  void cleanHisto();
     
-private:
   DetectorConstruction*   detector;
   PhysicsList*            physics;
@@ -84,8 +79,7 @@
   G4double                edeptot, eniel;
   G4int                   nPrimarySteps;
-           
-  AIDA::IAnalysisFactory* af;  
-  AIDA::ITree*            tree;
-  AIDA::IHistogram1D*     histo[1];        
+  G4int                   nRange;
+
+  Histo*                  histo;
 };
 
Index: trunk/examples/extended/electromagnetic/TestEm7/include/StepMax.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/StepMax.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/StepMax.hh	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: StepMax.hh,v 1.3 2006/06/29 16:58:01 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/StepMaxMessenger.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/StepMaxMessenger.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/StepMaxMessenger.hh	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: StepMaxMessenger.hh,v 1.2 2006/06/29 16:58:03 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/SteppingAction.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/SteppingAction.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/SteppingAction.hh	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: SteppingAction.hh,v 1.2 2006/06/29 16:58:05 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/SteppingVerbose.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/SteppingVerbose.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/SteppingVerbose.hh	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: SteppingVerbose.hh,v 1.2 2006/06/29 16:58:07 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/TrackingAction.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/TrackingAction.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/TrackingAction.hh	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: TrackingAction.hh,v 1.2 2006/06/29 16:58:09 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm7/include/c2_function.hh
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/c2_function.hh	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/c2_function.hh	(revision 1230)
@@ -1,2 +1,28 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+//
 /**
  *  \file 
@@ -7,13 +33,29 @@
  *  \author Copyright 2005 __Vanderbilt University__. All rights reserved.
  *
- * 	\version c2_function.hh,v 1.53 2007/11/12 13:58:57 marcus Exp
+ * 	\version c2_function.hh,v 1.238 2008/05/22 12:45:19 marcus Exp
+ *  \see \ref c2_factory "Factory Functions" for information on constructing things in here  
  */
 
-#ifndef __has_C2Functions_c2_h
-#define __has_C2Functions_c2_h 1
+#ifndef __has_c2_function_hh
+#define __has_c2_function_hh 1
+
+// MSVC does not automatically define numerical constants such as M_PI without this.
+// this came from the msdn website, so it should be right...
+#ifdef _MSC_VER
+#define _USE_MATH_DEFINES
+#define c2_isnan _isnan
+#define c2_isfinite _finite
+#else
+#define c2_isnan std::isnan
+#define c2_isfinite std::isfinite
+#endif
 
 #include <cmath>
 #include <vector>
+#include <utility>
 #include <string>
+#include <stdexcept>
+#include <typeinfo>
+#include <sstream>
 
 /// \brief the exception class for c2_function operations.
@@ -32,16 +74,49 @@
 
 // put these forward references here, and with a bogus typename to make swig happy.
-template <typename float_type> class c2_composed_function;
-template <typename float_type> class c2_sum;
-template <typename float_type> class c2_diff;
-template <typename float_type> class c2_product;
-template <typename float_type> class c2_ratio;
+template <typename float_type> class c2_composed_function_p;
+template <typename float_type> class c2_sum_p;
+template <typename float_type> class c2_diff_p;
+template <typename float_type> class c2_product_p;
+template <typename float_type> class c2_ratio_p;
+template <typename float_type> class c2_piecewise_function_p;
+template <typename float_type> class c2_quadratic_p;
+template <typename float_type> class c2_ptr;
+/**
+	\defgroup abstract_classes Abstract Classes
+	\defgroup arithmetic_functions Arithmetic Functions
+	\defgroup math_functions Mathemetical Functions
+	\defgroup parametric_functions Parametric Families of Functions
+	\defgroup interpolators Interpolating Functions
+	\defgroup containers Functions which are containers for, or functions of, other functions
+	\defgroup factories Factory classes which reduce silly template typing
+	\defgroup transforms Classes which provide coordinate system transformations, wih derivatives
+*/
+
+/// \brief structure used to hold evaluated function data at a point.  
+///
+/// Contains all the information for the function at one point. 
+template <typename float_type> class c2_fblock 
+{	
+public:
+	/// \brief the abscissa
+	float_type x;
+	/// \brief the value of the function at \a x
+	float_type y;
+	/// \brief the derivative at \a x
+	float_type yp;
+	/// \brief the second derivative at \a x
+	float_type ypp; 
+	/// flag, filled in by c2_function::fill_fblock(), indicating the derivative is NaN of Inf
+	bool ypbad;
+	/// flag, filled in by c2_function::fill_fblock(), indicating the second derivative is NaN of Inf
+	bool yppbad; 
+};
 
 /**
  \brief the parent class for all c2_functions.
-
+ \ingroup abstract_classes
   c2_functions know their value, first, and second derivative at almost every point.
   They can be efficiently combined with binary operators, via c2_binary_function, 
-  composed via c2_composed_function,
+  composed via c2_composed_function_,
   have their roots found via find_root(),
   and be adaptively integrated via partial_integrals() or integral().
@@ -53,15 +128,7 @@
     log_log_interpolating_function, and arrhenius_interpolating_function,
     as well as the template functions
-    inverse_integrated_density().
+    inverse_integrated_density_function().
  
- \warning
- The composite flavors of c2_functions (c2_sum, c2_composed_function, c2_binary_function, e.g.) make no effort to manage 
- deletion of their component functions.
- These are just container classes, and the user (along with normal automatic variable semantics) 
- is responsible for the lifetime of components.
- Inappropriate attention to this can cause massive memory leaks.  
- However, in most cases these do exactly what is intended.
- The classes will be left this way since the only other option is to use copy constructors on everything, 
- which would make this all very slow.
+ For a discussion of memory management, see \ref memory_management
  
  */
@@ -71,5 +138,5 @@
     /// \return the CVS Id string
 	const std::string cvs_header_vers() const { return 
-		"c2_function.hh,v 1.53 2007/11/12 13:58:57 marcus Exp";
+		"c2_function.hh,v 1.238 2008/05/22 12:45:19 marcus Exp";
 	}
 	
@@ -80,5 +147,14 @@
 public:
     /// \brief destructor
-	virtual ~c2_function() { if(sampling_grid && !no_overwrite_grid) delete sampling_grid; }
+	virtual ~c2_function() { 
+		if(sampling_grid && !no_overwrite_grid) delete sampling_grid; 	
+		if(root_info) delete root_info; 
+		if(owner_count) {
+			std::ostringstream outstr;
+			outstr << "attempt to delete an object with non-zero ownership in class ";
+			outstr << typeid(*this).name() << std::endl;
+			throw c2_exception(outstr.str().c_str());
+		}
+	}
 	
 	/// \brief get the value and derivatives. 
@@ -97,10 +173,4 @@
 	inline float_type operator () (float_type x) const throw(c2_exception) 
         { return value_with_derivatives(x, (float_type *)0, (float_type *)0); } 
-
-    /// \brief compose this function outside another.
-    /// \param inner the inner function
-    /// \return the composed function
-	c2_composed_function<float_type> & operator ()(const c2_function<float_type> &inner) const 
-		{ return *new c2_composed_function<float_type>((*this), inner); }
 
 	/// \brief get the value and derivatives. 
@@ -130,4 +200,5 @@
     /// \param[out] final_yprime2 If pointer is not zero, return second derivative of function at root
     /// \return the position of the root.
+	/// \see \ref rootfinder_subsec "Root finding sample" 
 	float_type find_root(float_type lower_bracket, float_type upper_bracket, float_type start, 
         float_type value, int *error=0, 
@@ -151,7 +222,8 @@
     /// with no error checking.
     /// \param extrapolate if true, use simple Richardson extrapolation on the final 2 steps to reduce the error.
-    /// \return sum of partial integrals, whcih is the definite integral from the first value in \a xgrid to the last.
+    /// \return sum of partial integrals, which is the definite integral from the first value in \a xgrid to the last.
 	float_type partial_integrals(std::vector<float_type> xgrid, std::vector<float_type> *partials = 0,
-          float_type abs_tol=1e-12, float_type rel_tol=1e-12, int derivs=2, bool adapt=true, bool extrapolate=true) const;
+          float_type abs_tol=1e-12, float_type rel_tol=1e-12, int derivs=2, bool adapt=true, bool extrapolate=true) 
+		const throw(c2_exception);
 	
     /// \brief a fully-automated integrator which uses the information provided by the get_sampling_grid() function
@@ -164,5 +236,5 @@
 	/// \param xmax upper bound of the domain for integration
     /// \param partials if non-NULL, a vector in which to receive the partial integrals.
-	/// It will automatically be sized apprpropriately, if provided, to contain \a n - 1 elements where \a n is the length of \a xgrid  
+	/// It will automatically be sized appropriately, if provided, to contain \a n - 1 elements where \a n is the length of \a xgrid  
     /// \param abs_tol the absolute error bound for each segment
     /// \param rel_tol the fractional error bound for each segment.  
@@ -173,8 +245,50 @@
     /// with no error checking.
     /// \param extrapolate if true, use simple Richardson extrapolation on the final 2 steps to reduce the error.
-    /// \return sum of partial integrals, whcih is the definite integral from the first value in \a xgrid to the last.
+    /// \return sum of partial integrals, which is the definite integral from the first value in \a xgrid to the last.
 	float_type integral(float_type xmin, float_type xmax, std::vector<float_type> *partials = 0,
-             float_type abs_tol=1e-12, float_type rel_tol=1e-12, int derivs=2, bool adapt=true, bool extrapolate=true) const;
-
+             float_type abs_tol=1e-12, float_type rel_tol=1e-12, int derivs=2, bool adapt=true, bool extrapolate=true) 
+		const throw(c2_exception);
+
+	/// \brief create a c2_piecewise_function_p from c2_connector_function_p segments which 
+	/// is a representation of the parent function to the specified accuracy, but maybe much cheaper to evaluate
+	///
+	/// This method has three modes, depending on the \a derivs flag. 
+	///
+	/// If \a derivs is 2,
+	/// it computes a c2_piecewise_function_p representation of its parent function, which may be a much faster 
+	/// function to use in codes if the parent function is expensive.  If \a xvals and \a yvals are non-null,
+	/// it will also fill them in with the function values at each grid point the adaptive algorithm chooses.
+	///
+	/// If \a derivs is 1, this does not create the connectors, 
+	/// and returns an null pointer, but will fill in the \a xvals and \a yvals 
+	/// vectors with values of the function at points such that the linear interpolation error between the points
+	/// is bounded by the tolerance values given.  Because it uses derivative information from the function to manage the 
+	/// error control, it is almost completely free of issues with missing periods of oscillatory functions,
+	/// even with no information provided in the sampling grid.
+	/// This is typically useful for sampling a function for plotting.
+	///
+	/// If \a derivs is 0, this does something very like what it does if \a derivs = 1, but without derivatives.  
+	/// Instead, to compute the intermediate value of the function for error control, it just uses 
+	/// 3-point parabolic interpolation.  This is useful amost exclusively for converting a non-c2_function,
+	/// with no derivatives, but wrapped in a c2_classic_function wrapper, into a table of values to seed an interpolating_function_p.
+	/// Note, however, that without derivatives, this is very susceptible to missing periods of oscillatory 
+	/// functions, so it is important to set a sampling grid which isn't too much coarser than the typical oscillations.
+	///
+	/// \note the \a sampling_grid of the returned function matches the \a sampling_grid of its parent.
+	/// \see \ref sample_function_for_plotting "Adaptive Sampling Examples"
+    /// \param xmin lower bound of the domain for sampling
+	/// \param xmax upper bound of the domain for sampling
+    /// \param abs_tol the absolute error bound for each segment
+    /// \param rel_tol the fractional error bound for each segment. 
+	/// \param derivs if 0 or 1, return a useless function, but fill in the \a xvals and \a yvals vectors (if non-null).
+	/// Also, if 0 or 1, tolerances refer to linear interpolation, not high-order interpolation. 
+	/// If 2, return a full piecewise collection of c2_connector_function_p segments.  See discussion above.
+	/// \param [in,out] xvals vector of abscissas at which the function was actually sampled (if non-null)
+	/// \param [in,out] yvals vector of function values corresponding to \a xvals (if non-null)
+	/// \return a new, sampled representation, if \a derivs is 2.  A null pointer if \a derivs is 0 or 1.
+	c2_piecewise_function_p<float_type> *adaptively_sample(float_type xmin, float_type xmax,
+		 float_type abs_tol=1e-12, float_type rel_tol=1e-12,
+		 int derivs=2, std::vector<float_type> *xvals=0, std::vector<float_type> *yvals=0) const throw(c2_exception);
+	
 	/// \brief return the lower bound of the domain for this function as set by set_domain()
 	inline float_type xmin() const { return fXMin; }
@@ -186,8 +300,7 @@
 	/// \brief this is a counter owned by the function but which can be used to monitor efficiency of algorithms.
 	///
-	/// It is not maintained automatically in general!  
-	/// The adaptive integrator and root finder do clear it at the start and update it for performance checking.
+	/// It is not maintained automatically in general!  The root finder, integrator, and sampler do increment it.
 	/// \return number of evaluations logged since last reset.
-	volatile int get_evaluations() const { return evaluations; }
+	volatile size_t get_evaluations() const { return evaluations; }
 	/// \brief reset the counter
 	void reset_evaluations()  const { evaluations=0; } // evaluations are 'invisible' to constant
@@ -200,5 +313,5 @@
 	/// \param message an informative string to include in an exception if this throws c2_exception
 	/// \return true if in decreasing order, false if increasing 
-	bool check_monotonicity(const std::vector<float_type> &data, const char message[]) throw(c2_exception);
+	bool check_monotonicity(const std::vector<float_type> &data, const char message[]) const throw(c2_exception);
 	
 	/// \brief establish a grid of 'interesting' points on the function.
@@ -212,4 +325,8 @@
 	virtual void set_sampling_grid(const std::vector<float_type> &grid) throw(c2_exception); 
 	
+	/// \brief get the sampling grid, which may be a null pointer
+	/// \return pointer to the sampling grid
+	std::vector<float_type> *get_sampling_grid_pointer() const { return sampling_grid; } 
+	
     /// \brief return the grid of 'interesting' points along this function which lie in the region requested
 	///
@@ -217,17 +334,15 @@
 	/// \param xmin the lower bound for which the function is to be sampled
 	/// \param xmax the upper bound for which the function is to be sampled
-	/// \return a new vector containing the samplng grid.
-	virtual std::vector<float_type> &get_sampling_grid(float_type xmin, float_type xmax) const ;
+	/// \param [in,out] grid filled vector containing the samplng grid.
+	virtual void get_sampling_grid(float_type xmin, float_type xmax, std::vector<float_type> &grid) const ;
 	
 	/// \brief clean up endpoints on a grid of points
-	///
 	/// \param[in,out] result the sampling grid with excessively closely space endpoints removed.
 	/// The grid is modified in place.
 	void preen_sampling_grid(std::vector<float_type> *result) const;		
 	/// \brief refine a grid by splitting each interval into more intervals
-	/// \param grid the grid to refine
+	/// \param [in,out] grid the grid to refine in place
 	/// \param refinement the number of new steps for each old step
-	/// \return a new sampling grid with more points.
-	std::vector<float_type> & refine_sampling_grid(const std::vector<float_type> &grid, size_t refinement) const;		
+	void refine_sampling_grid(std::vector<float_type> &grid, size_t refinement) const;		
 	
 	/// \brief create a new c2_function from this one which is normalized on the interval 
@@ -236,5 +351,5 @@
 	/// \param norm the desired integral for the function over the region
 	/// \return a new c2_function with the desired \a norm.
-	c2_function<float_type> &normalized_function(float_type xmin, float_type xmax, float_type norm=1.0);
+	c2_function<float_type> &normalized_function(float_type xmin, float_type xmax, float_type norm=1.0) const throw(c2_exception);
 	/// \brief create a new c2_function from this one which is square-normalized on the interval 
     /// \param xmin lower bound of the domain for integration
@@ -242,5 +357,5 @@
 	/// \param norm the desired integral for the function over the region
 	/// \return a new c2_function with the desired \a norm.
-	c2_function<float_type> &square_normalized_function(float_type xmin, float_type xmax, float_type norm=1.0);
+	c2_function<float_type> &square_normalized_function(float_type xmin, float_type xmax, float_type norm=1.0) const throw(c2_exception);
 	/// \brief create a new c2_function from this one which is square-normalized with the provided \a weight on the interval 
     /// \param xmin lower bound of the domain for integration
@@ -249,38 +364,71 @@
 	/// \param norm the desired integral for the function over the region
 	/// \return a new c2_function with the desired \a norm.
-	c2_function<float_type> &square_normalized_function(float_type xmin, float_type xmax, const c2_function<float_type> &weight, float_type norm=1.0);
-
-	/// \brief factory function to create a c2_sum from an regular algebraic expression.
-	/// \note
-	/// be very wary of ownership issues if this is used in a complex expression.
-	c2_sum<float_type> &operator + (const c2_function<float_type> &rhs)  { return *new c2_sum<float_type>(*this, rhs); }
-	/// \brief factory function to create a c2_diff from an regular algebraic expression.
-	/// \note
-	/// be very wary of ownership issues if this is used in a complex expression.
-	c2_diff<float_type> &operator - (const c2_function<float_type> &rhs)  { return *new c2_diff<float_type>(*this, rhs); }
-	/// \brief factory function to create a c2_product from an regular algebraic expression.
-	/// \note
-	/// be very wary of ownership issues if this is used in a complex expression.
-	c2_product<float_type> &operator * (const c2_function<float_type> &rhs)  { return *new c2_product<float_type>(*this, rhs); }
-	/// \brief factory function to create a c2_ratio from an regular algebraic expression.
-	/// \note
-	/// be very wary of ownership issues if this is used in a complex expression.
-	c2_ratio<float_type> &operator / (const c2_function<float_type> &rhs)  { return *new c2_ratio<float_type>(*this, rhs); }
-	
-
-	
-	std::vector<float_type> *sampling_grid;
-	bool no_overwrite_grid;
-	    
+	c2_function<float_type> &square_normalized_function(
+		float_type xmin, float_type xmax, const c2_function<float_type> &weight, float_type norm=1.0)
+		const throw(c2_exception);
+
+	/// \brief factory function to create a c2_sum_p from a regular algebraic expression.
+	/// \param rhs the right-hand term of the sum
+	/// \return a new c2_function
+	c2_sum_p<float_type> &operator + (const c2_function<float_type> &rhs)  const 
+		{ return *new c2_sum_p<float_type>(*this, rhs); }
+	/// \brief factory function to create a c2_diff_p from a regular algebraic expression.
+	/// \param rhs the right-hand term of the difference
+	/// \return a new c2_function
+	c2_diff_p<float_type> &operator - (const c2_function<float_type> &rhs)  const 
+		{ return *new c2_diff_p<float_type>(*this, rhs); }
+	/// \brief factory function to create a c2_product_p from a regular algebraic expression.
+	/// \param rhs the right-hand term of the product
+	/// \return a new c2_function
+	c2_product_p<float_type> &operator * (const c2_function<float_type> &rhs) const 
+		{ return *new c2_product_p<float_type>(*this, rhs); }
+	/// \brief factory function to create a c2_ratio_p from a regular algebraic expression.
+	/// \param rhs the right-hand term of the ratio (the denominator)
+	/// \return a new c2_function
+	c2_ratio_p<float_type> &operator / (const c2_function<float_type> &rhs)  const 
+		{ return *new c2_ratio_p<float_type>(*this, rhs); }
+    /// \brief compose this function outside another.
+    /// \param inner the inner function
+    /// \return the composed function
+	/// \anchor compose_operator
+	c2_composed_function_p<float_type> & operator ()(const c2_function<float_type> &inner) const 
+		{ return *new c2_composed_function_p<float_type>((*this), inner); }
+
+	/// \brief Find out where a calculation ran into trouble, if it got a nan.
+	/// If the most recent computation did not return a nan, this is undefined.
+	/// \return \a x value of point at which something went wrong, if integrator (or otherwise) returned a nan.
+	float_type get_trouble_point() const { return bad_x_point; }
+	
+	/// \brief increment our reference count.  Destruction is only legal if the count is zero.
+	void claim_ownership() const { owner_count++; }
+	/// \brief decrement our reference count. Do not destroy at zero.
+	/// \return final owner count, to check whether object should disappear.
+	size_t release_ownership_for_return() const throw(c2_exception) { 
+		if(!owner_count) {
+			std::ostringstream outstr;
+			outstr << "attempt to release ownership of an unowned function in class ";
+			outstr << typeid(*this).name() << std::endl;
+			throw c2_exception(outstr.str().c_str());
+		}
+		owner_count--; 
+		return owner_count;
+	}
+	/// \brief decrement our reference count. If the count reaches zero, destroy ourself.
+	void release_ownership() const throw(c2_exception) { 
+		if(!release_ownership_for_return()) delete this;
+	}
+	/// \brief get the reference count, mostly for debugging
+	/// \return the count
+	size_t count_owners() const { return owner_count; }
+
 protected:
 	c2_function(const c2_function<float_type>  &src) : sampling_grid(0),
 		no_overwrite_grid(false),
-        fXMin(src.fXMin), fXMax(src.fXMax), rootInitialized(false)
+        fXMin(src.fXMin), fXMax(src.fXMax), root_info(0), owner_count(0)
         {} // copy constructor only copies domain, and is only for internal use
 	c2_function() : 
 			sampling_grid(0), no_overwrite_grid(0), 
         fXMin(-std::numeric_limits<float_type>::max()), 
-        fXMax(std::numeric_limits<float_type>::max()),
-        rootInitialized(false)
+        fXMax(std::numeric_limits<float_type>::max()), root_info(0), owner_count(0)
 		{} // prevent accidental naked construction (impossible any since this has pure virtual methods)
 	
@@ -293,20 +441,61 @@
 		}
 	
+	std::vector<float_type> * sampling_grid;
+	bool no_overwrite_grid;
+	
 	float_type fXMin, fXMax;
-	mutable int evaluations;
-	
+	mutable size_t evaluations;
+	/// \brief this point may be used to record where a calculation ran into trouble
+	mutable float_type bad_x_point;
+public: 
+	/// \brief fill in a c2_fblock<float_type>... a shortcut for the integrator & sampler
+	/// \param [in,out] fb the block to fill in with information
+	inline void fill_fblock(c2_fblock<float_type> &fb) const throw(c2_exception) 
+	{
+		fb.y=value_with_derivatives(fb.x, &fb.yp, &fb.ypp);
+		fb.ypbad=c2_isnan(fb.yp) || !c2_isfinite(fb.yp);
+		fb.yppbad=c2_isnan(fb.ypp) || !c2_isfinite(fb.ypp);
+		increment_evaluations();
+	}
+
 private:
-	/// \brief structure used for recursion in adaptive integrator.  
-	///
-	/// Contains all the information for the function at one point. 
-	struct c2_integrate_fblock {	float_type x, y, yp, ypp; };
-	/// \brief structure used to pass information recursively.
+	/// \brief the data element for the internal recursion stack for the sampler and integrator
+	struct recur_item { 
+		c2_fblock<float_type> f1; size_t depth; 
+		float_type previous_estimate, abs_tol, step_sum; 
+		bool done; 
+		size_t f0index, f2index;
+	};
+	
+
+	/// \brief structure used to pass information recursively in integrator.
 	///
 	/// the \a abs_tol is scaled by a factor of two at each division.  
 	/// Everything else is just passed down.
-	struct c2_integrate_recur { struct c2_integrate_fblock *f0, *f1, *f2;
-		float_type abs_tol, rel_tol, *lr, eps_scale, extrap_coef, extrap2;
-		int depth, derivs;
-		bool adapt, extrapolate;
+	struct c2_integrate_recur { 
+		c2_fblock<float_type> *f0, *f1;
+		float_type abs_tol, rel_tol, eps_scale, extrap_coef, extrap2, dx_tolerance, abs_tol_min;
+		std::vector< recur_item > *rb_stack;
+		int  derivs;
+		bool adapt, extrapolate, inited;
+	};
+
+	/// \brief structure used to pass information recursively in sampler.
+	///
+	struct c2_sample_recur { 
+		c2_fblock<float_type> *f0, *f1;
+		float_type abs_tol, rel_tol, dx_tolerance, abs_tol_min;
+		int derivs;
+		c2_piecewise_function_p<float_type> *out;
+		std::vector<float_type> *xvals, *yvals;
+		std::vector< recur_item > *rb_stack;
+		bool inited;
+	};
+
+	/// \brief structure used to hold root bracketing information
+	///
+	struct c2_root_info {
+		c2_fblock<float_type> lower, upper;
+		bool inited;
 	};
 	
@@ -316,17 +505,265 @@
     /// to allow very efficient recursion.
     /// \param rb a pointer to the recur struct.
-	float_type integrate_step(struct c2_integrate_recur &rb) const;
+	float_type integrate_step(struct c2_integrate_recur &rb) const throw(c2_exception);
     
-    /// these carry a memory of the last root bracketing,
+    /// \brief Carry out the recursive subdivision for sampling.
+    ///
+    /// This passes information recursively through the \a recur block pointer
+    /// to allow very efficient recursion.
+    /// \param rb a pointer to the recur struct.
+	void sample_step(struct c2_sample_recur &rb) const throw(c2_exception);
+
+    /// this carry a memory of the last root bracketing,
     /// to avoid the necessity of evaluating the function on the brackets every time
     /// if the brackets have not been changed. 
-    mutable float_type lastRootLowerX, lastRootUpperX, lastRootLowerY, lastRootUpperY;
-    mutable int rootInitialized;
-	
-};
-
+	/// it is declared as a pointer, since many c2_functions may never need one allocated
+    mutable struct  c2_root_info *root_info;
+	
+	mutable size_t owner_count;
+};
+
+/// \brief a container into which any conventional c-style function can be dropped,
+/// to create a degenerate c2_function without derivatives. 
+/// Mostly useful for sampling into interpolating functions.
+/// construct a reference to this with c2_classic_function()
+/// \ingroup containers
+/// The factory function c2_factory::classic_function() creates *new c2_classic_function_p()
+template <typename float_type=double> class c2_classic_function_p : public c2_function<float_type> {
+public:
+	/// \brief construct the container 
+	/// \param c_func a pointer to a conventional c-style function
+	c2_classic_function_p(const float_type (*c_func)(float_type)) : c2_function<float_type>(), func(c_func)  {}
+
+	/// \copydoc c2_function::value_with_derivatives
+	/// Uses the internal function pointer set by set_function().
+	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) 
+	{
+		if(!func) throw c2_exception("c2_classic_function called with null function");
+		if(yprime) *yprime=0;
+		if(yprime2) *yprime2=0;
+		return func(x); 
+	}
+	~c2_classic_function_p() { }
+	
+protected:
+	/// \brief pointer to our function
+	const float_type (*func)(float_type);
+};
+
+/// \brief create a container for a c2_function which handles the reference counting.
+/// \ingroup containers
+/// It is useful as a smart container to hold a c2_function and keep the reference count correct.  
+/// The recommended way for a class to store a c2_function which is handed in from the outside 
+/// is for it to have a c2_ptr member into which the passed-in function is stored.
+/// This way, when the class instance is deleted, it will automatically dereference any function
+/// which it was handed.
+///
+/// This class contains a copy constructor and operator=, to make it fairly easy to make 
+/// a std::vector of these objects, and have it work as expected.
+template <typename float_type> class c2_const_ptr {
+public:
+	/// \brief construct the container with no function
+	c2_const_ptr() : func(0)  {}
+	/// \brief construct the container with a pre-defined function
+	/// \param f the function to store
+	c2_const_ptr(const c2_function<float_type> &f) : func(0) 
+		{ set_function(&f); }
+	/// \brief copy constructor
+	/// \param src the container to copy
+	c2_const_ptr(const c2_const_ptr<float_type> &src) : func(0)
+		{ set_function(src.get_ptr()); }
+	/// \brief fill the container with a new function, or clear it with a null pointer
+	/// \param f the function to store, releasing any previously held function
+	void set_function(const c2_function<float_type> *f) 
+	{ 
+		if(func) func->release_ownership();
+		func=f; 
+		if(func) func->claim_ownership();
+	}
+	
+	/// \brief fill the container from another container
+	/// \param f the container to copy
+	void operator =(const c2_const_ptr<float_type> &f) 
+		{ set_function(f.get_ptr()); }
+	/// \brief fill the container with a function
+	/// \param f the function
+	void operator =(const c2_function<float_type> &f) 
+		{ set_function(&f); }
+	/// \brief release the function without destroying it, so it can be returned from a function
+	///
+	/// This is usually the very last line of a function before the return statement, so that 
+	/// any exceptions that happen during execution of the function will cause proper cleanup.
+	/// Once the function has been released from its container this way, it is an orhpaned object
+	/// until the caller claims it, so it could get lost if an exception happens.
+	void release_for_return() throw(c2_exception)
+		{	
+			if(func) func->release_ownership_for_return();
+			func=0;
+		}
+	/// \brief clear the function
+	///
+	/// Any attempt to use this c2_plugin_function_p throws an exception if the saved function is cleared.
+	void unset_function(void) { set_function(0);  }
+	/// \brief destructor
+	~c2_const_ptr() { set_function(0); }
+	
+	/// \brief get a reference to our owned function
+	inline const c2_function<float_type> &get() const throw(c2_exception) 
+	{ 
+		if(!func) throw c2_exception("c2_ptr accessed uninitialized");
+		return *func;
+	}
+	/// \brief get an unchecked pointer to our owned function
+	inline const c2_function<float_type> *get_ptr() const { return func; }
+	/// \brief get a checked pointer to our owned function
+	inline const c2_function<float_type> *operator -> () const 
+		{ return &get(); }
+	/// \brief check if we have a valid function
+	bool valid() const { return func != 0; }
+	
+	/// \brief type coercion operator which lets us use a pointer as if it were a const c2_function
+	operator const c2_function<float_type>& () const { return this->get(); }
+	
+	/// \brief convenience operator to make us look like a function
+	/// \param x the value at which to evaluate the contained function
+	/// \return the evaluated function
+	/// \note If you using this repeatedly, do const c2_function<float_type> &func=ptr;
+	/// and use func(x).  Calling this operator wastes some time, since it checks the validity of the
+	/// pointer every time.
+	float_type operator()(float_type x) const throw(c2_exception) { return get()(x); }
+	/// \brief convenience operator to make us look like a function
+	/// \param x the value at which to evaluate the contained function
+	/// \param yprime the derivative
+	/// \param yprime2 the second derivative
+	/// \return the evaluated function
+	/// \note If you using this repeatedly, do const c2_function<float_type> &func=ptr;
+	/// and use func(x).  Calling this operator wastes some time, since it checks the validity of the
+	/// pointer every time.
+	float_type operator()(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) 
+		{ return get().value_with_derivatives(x, yprime, yprime2); }
+	/// \brief factory function to create a c2_sum_p from a regular algebraic expression.
+	/// \param rhs the right-hand term of the sum
+	/// \return a new c2_function
+	c2_sum_p<float_type> &operator + (const c2_function<float_type> &rhs)  const throw(c2_exception)
+		{ return *new c2_sum_p<float_type>(get(), rhs); }
+	/// \brief factory function to create a c2_diff_p from a regular algebraic expression.
+	/// \param rhs the right-hand term of the difference
+	/// \return a new c2_function
+	c2_diff_p<float_type> &operator - (const c2_function<float_type> &rhs)  const throw(c2_exception)
+		{ return *new c2_diff_p<float_type>(get(), rhs); }
+	/// \brief factory function to create a c2_product_p from a regular algebraic expression.
+	/// \param rhs the right-hand term of the product
+	/// \return a new c2_function
+	c2_product_p<float_type> &operator * (const c2_function<float_type> &rhs) const throw(c2_exception)
+		{ return *new c2_product_p<float_type>(get(), rhs); }
+	/// \brief factory function to create a c2_ratio_p from a regular algebraic expression.
+	/// \param rhs the right-hand term of the ratio (the denominator)
+	/// \return a new c2_function
+	c2_ratio_p<float_type> &operator / (const c2_function<float_type> &rhs)  const throw(c2_exception)
+		{ return *new c2_ratio_p<float_type>(get(), rhs); }
+    /// \brief compose this function outside another.
+    /// \param inner the inner function
+    /// \return the composed function
+	c2_composed_function_p<float_type> & operator ()(const c2_function<float_type> &inner) const throw(c2_exception)
+		{ return *new c2_composed_function_p<float_type>(get(), inner); }
+	
+protected:
+	const c2_function<float_type> * func;
+};
+
+/// \brief create a container for a c2_function which handles the reference counting.
+/// \ingroup containers
+///
+/// \see  c2_const_ptr and \ref memory_management "Use of c2_ptr for memory management"
+
+template <typename float_type> class c2_ptr : public c2_const_ptr<float_type >
+{
+public:
+	/// \brief construct the container with no function
+	c2_ptr() : c2_const_ptr<float_type>()  {}
+	/// \brief construct the container with a pre-defined function
+	/// \param f the function to store
+	c2_ptr(c2_function<float_type> &f) : 
+		c2_const_ptr<float_type>() { set_function(&f); } 
+	/// \brief copy constructor
+	/// \param src the container to copy
+	c2_ptr(const c2_ptr<float_type> &src) : 
+		c2_const_ptr<float_type>() { set_function(src.get_ptr()); }
+	/// \brief get a checked pointer to our owned function
+	inline c2_function<float_type> &get() const throw(c2_exception) 
+		{ return *const_cast<c2_function<float_type>*>(&c2_const_ptr<float_type>::get()); }
+	/// \brief get an unchecked pointer to our owned function
+	inline c2_function<float_type> *get_ptr() const 
+		{ return const_cast<c2_function<float_type>*>(this->func); }
+	/// \brief get a checked pointer to our owned function
+	inline c2_function<float_type> *operator -> () const 
+		{ return &get(); }
+	/// \brief fill the container from another container
+	/// \param f the container to copy
+	void operator =(const c2_ptr<float_type> &f) 
+		{ set_function(f.get_ptr()); }
+	/// \brief fill the container with a function
+	/// \param f the function
+	void operator =(c2_function<float_type> &f) 
+		{ set_function(&f); }
+private:
+	/// \brief hidden non-const-safe version of operator=
+	void operator =(const c2_const_ptr<float_type> &f) { }
+	/// \brief hidden non-const-safe version of operator=
+	void operator =(const c2_function<float_type> &f) { }
+};
+
+/// \brief create a non-generic container for a c2_function which handles the reference counting.
+/// \ingroup containers
+///
+/// \see  c2_const_ptr and \ref memory_management "Use of c2_ptr for memory management"
+///
+/// \note Overuse of this class will generate massive bloat.  Use c2_ptr and c2_const_ptr if you don't _really_ need specific pointer types.
+/// \see  \ref memory_management "Use of c2_ptr for memory management"
+/*
+template <typename float_type, template <typename> class c2_class > class c2_typed_ptr : public c2_const_ptr<float_type> {
+public:
+	/// \brief construct the container with no function
+	c2_typed_ptr() : c2_ptr<float_type>()  {}
+	/// \brief construct the container with a pre-defined function
+	/// \param f the function to store
+	c2_typed_ptr(c2_class<float_type> &f)
+		: c2_const_ptr<float_type>() { this->set_function(&f); }
+	/// \brief copy constructor
+	/// \param src the container to copy
+	c2_typed_ptr(const c2_typed_ptr<float_type, c2_class> &src) 
+		: c2_const_ptr<float_type>() { this->set_function(src.get_ptr()); }
+	
+	/// \brief get a reference to our owned function
+	inline c2_class<float_type> &get() const throw(c2_exception) 
+		{ 
+			return *static_cast<c2_class<float_type> *>(const_cast<c2_function<float_type>*>(&c2_const_ptr<float_type>::get()));
+		}
+	/// \brief get a checked pointer to our owned function
+	inline c2_class<float_type> *operator -> () const 
+		{ return &get(); }
+	/// \brief get an unchecked pointer to our owned function
+	inline c2_class<float_type> *get_ptr() const 
+		{ return static_cast<c2_class<float_type> *>(const_cast<c2_function<float_type>*>(this->func)); }
+	/// \brief type coercion operator which lets us use a pointer as if it were a c2_function
+	operator c2_class<float_type>&() const { return get(); }
+	/// \brief fill the container from another container
+	/// \param f the container to copy
+	void operator =(const c2_typed_ptr<float_type, c2_class> &f) 
+		{ set_function(f.get_ptr()); }
+	/// \brief fill the container with a function
+	/// \param f the function
+	void operator =(c2_class<float_type> &f) 
+		{ set_function(&f); }
+private:
+	/// \brief hidden downcasting version of operator=
+	void operator =(const c2_const_ptr<float_type> &f) { }
+	/// \brief hidden downcasting version of operator=. Use an explicit dynamic_cast<c2_class<float_type>&>(f) if you need to try this.
+	void operator =(const c2_function<float_type> &f) { }
+};
+*/
 /// \brief a container into which any other c2_function can be dropped, to allow expressions
 /// with replacable components.  
-///
+/// \ingroup containers
 ///It is useful for plugging different InterpolatingFunctions into a c2_function expression.
 ///It saves a lot of effort in other places with casting away const declarations.
@@ -334,20 +771,26 @@
 /// It is also useful as a wrapper for a function if it is necessary to have a copy of a function
 /// which has a different domain or sampling grid than the parent function.  This can be
-/// be used, for example, to patch badly-behaved functions with c2_piecewise_function by 
+/// be used, for example, to patch badly-behaved functions with c2_piecewise_function_p by 
 /// taking the parent function, creating two plugins of it with domains on each side of the 
 /// nasty bit, and then inserting a nice function in the hole.
-template <typename float_type=double> class c2_plugin_function : public c2_function<float_type> {
+///
+/// This can also be used as a fancier c2_ptr which allows direct evaluation
+/// instead of having to dereference the container first.
+///
+/// The factory function c2_factory::plugin_function() creates *new c2_plugin_function_p()
+template <typename float_type=double> class c2_plugin_function_p : 
+	public c2_function<float_type> {
 public:
 	/// \brief construct the container with no function
-	c2_plugin_function() : c2_function<float_type>(), func(0), owns(false)  {}
+	c2_plugin_function_p() : c2_function<float_type>(), func()  {}
 	/// \brief construct the container with a pre-defined function
-	c2_plugin_function(const c2_function<float_type> &f) : c2_function<float_type>(), func(0), owns(false) { set_function(f); }
-	/// \brief fill the container with a new function
-	void set_function(const c2_function<float_type> &f) 
-		{ 
-			if(owns && func) delete func;
-			func=&f; set_domain(f.xmin(), f.xmax()); 
-			set_sampling_grid_pointer(*f.sampling_grid); 
-			owns=false;
+	c2_plugin_function_p(c2_function<float_type> &f) : 
+		c2_function<float_type>(),func(f)  { }
+	/// \brief fill the container with a new function, or clear it with a null pointer
+	/// and copy our domain
+	void set_function(c2_function<float_type> *f) 
+		{
+			func.set_function(f);
+			if(f) set_domain(f->xmin(), f->xmax());
 		}
 	/// \copydoc c2_function::value_with_derivatives
@@ -355,100 +798,103 @@
 	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) 
 		{
-			if(!func) throw c2_exception("c2_plugin_function<float_type> called uninitialized");
-			return this->func->value_with_derivatives(x, yprime, yprime2); 
+			if(!func.valid()) throw c2_exception("c2_plugin_function_p called uninitialized");
+			return func->value_with_derivatives(x, yprime, yprime2); 
 		}
-	/// \brief clear the function
-	///
-	/// Any attempt to use this c2_plugin_function throws an exception if the saved function is cleared.
-	void unset_function(void) { if(owns && func) delete func; func=0; owns=false; }
 	/// \brief destructor
-	~c2_plugin_function() { if(owns && func) delete func; }
-	/// \brief tell us we should delete the function at destruction.  NOT sticky when function is reset
-	void set_ownership() { this->owns=true; }
-	
+	~c2_plugin_function_p() { }
+	
+	/// \brief clear our function
+	void unset_function() { func.unset_function(); }
+	
+	virtual void get_sampling_grid(float_type xmin, float_type xmax, std::vector<float_type> &grid) const 
+	{	
+		if(!func.valid()) throw c2_exception("c2_plugin_function_p called uninitialized");
+		if(this->sampling_grid) c2_function<float_type>::get_sampling_grid(xmin, xmax, grid);
+		else  func->get_sampling_grid(xmin, xmax, grid);
+	}
 protected:
-	const c2_function<float_type> *func;
-	bool owns;
-	
+	c2_ptr<float_type> func;
+};
+
+/// \brief a c2_plugin_function_p which promises not to fiddle with the plugged function.
+/// \ingroup containers
+///
+/// The factory function c2_factory::const_plugin_function() creates *new c2_const_plugin_function_p()
+template <typename float_type=double> class c2_const_plugin_function_p : public c2_plugin_function_p<float_type> {
+public:
+	/// \brief construct the container with no function
+	c2_const_plugin_function_p() : c2_plugin_function_p<float_type>()  {}
+	/// \brief construct the container with a pre-defined function
+	c2_const_plugin_function_p(const c2_function<float_type> &f) : 
+		c2_plugin_function_p<float_type>()  { set_function(&f); }
+	/// \brief fill the container with a new function, or clear it with a null pointer
+	void set_function(const c2_function<float_type> *f) 
+		{ c2_plugin_function_p<float_type>::set_function(const_cast<c2_function<float_type>*>(f)); }
+	/// \brief destructor
+	~c2_const_plugin_function_p() { }
+	
+	/// \brief get a const reference to our owned function, for direct access
+	const c2_function<float_type> &get() const throw(c2_exception) 
+		{ return this->func.get(); }
 };
 
 /// \brief Provides support for c2_function objects which are constructed from two other c2_function
 /// objects.  
-///
-/// It provides a very primitive ownership concept, so that the creator can tag a function
-/// as being owned by this function, so that when this function is deleted, the owned function will be deleted, too.
-/// This allows a piece of code to create various c2_function objects, combine them with binary operators,
-/// appropriately mark wich ones have no other possible owners, and return the final function with
-/// reasonable faith that everything will get cleaned up when the final function is deleted.  Note that
-/// none of this marking is automatic, to keep this class very lightweight.
+/// \ingroup abstract_classes
 template <typename float_type=double> class c2_binary_function : public c2_function<float_type> {
 public:
-	
-	
 	///  \brief function to manage the binary operation, used by c2_binary_function::value_with_derivatives() 
     /// 
-    /// Normally not used alone, but can be used to combine functions in other contexts. 
-	/// See interpolating_function::binary_operator() for an example.
-	/// \param left the function on the left of the binary operator or outside the composition
-	/// \param right the function to the right of the operator or inside the composition
-    /// \param[in] x the point at which to evaluate the function
-    /// \param[out] yprime the first derivative (if pointer is non-null)
-    /// \param[out] yprime2 the second derivative (if pointer is non-null)
-    /// \return the value of the function
-	virtual float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, 
-					   float_type x, float_type *yprime, float_type *yprime2) const =0;
 	
 	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception) 
     { 
-		return this->combine(this->Left, this->Right, x, yprime, yprime2);
+		if(stub) throw c2_exception("attempt to evaluate a c2_binary_function stub");
+		return this->combine(*Left.get_ptr(), *Right.get_ptr(), x, yprime, yprime2);
     }
 
-	/// \brief allow c2_binary_function to remember ownership of contained functions for automatic cleanup
+	/// \brief destructor releases ownership of member functions
 	///
-	/// upon destruction, this will cause disposal of the left member function
-	void set_left_ownership(void) { leftown=true; }
-	/// \brief allow c2_binary_function to remember ownership of contained functions for automatic cleanup
-	///
-	/// upon destruction, this will cause disposal of the right member function
-	void set_right_ownership(void) { rightown=true; }
-	/// \brief allow c2_binary_function to remember ownership of contained functions for automatic cleanup
-	///
-	/// upon destruction, this will cause disposal of both member functions
-	void set_ownership(void) { leftown=rightown=true; }
-	/// \brief destructor executes disposal of member functions if flagged
-	///
-	/// depends on judicious use of set_ownership(), set_right_ownership(), or set_left_ownership()
-	virtual ~c2_binary_function() {
-		if(leftown) delete &Left;
-		if(rightown) delete &Right;
-	}
-		
+	virtual ~c2_binary_function() { }
+
 protected:
 	/// \brief construct the binary function
+	/// \param combiner pointer to the function which actualy knows how to execute the binary
 	/// \param left the c2_function to be used in the left side of the binary relation
 	/// \param right the c2_function to be used in the right side of the binary relation
-	c2_binary_function( const c2_function<float_type> &left,  const c2_function<float_type> &right) : 
-	c2_function<float_type>(), 
-	Left(left), Right(right), leftown(false), rightown(false) 
+	c2_binary_function( 
+			float_type (*combiner)(const c2_function<float_type> &left, const c2_function<float_type> &right, 
+				   float_type x, float_type *yprime, float_type *yprime2),						
+			const c2_function<float_type> &left,  const c2_function<float_type> &right) : 
+			c2_function<float_type>(), combine(combiner), Left(left), Right(right), stub(false)
 	{ 
-			set_domain(
-					   (left.xmin() > right.xmin()) ? left.xmin() : right.xmin(), 
-					   (left.xmax() < right.xmax()) ? left.xmax() : right.xmax()
-					   );
+		set_domain(
+				   (left.xmin() > right.xmin()) ? left.xmin() : right.xmin(), 
+				   (left.xmax() < right.xmax()) ? left.xmax() : right.xmax()
+				   );
 	} 
 	
 	/// \brief construct a 'stub' c2_binary_function, which provides access to the combine() function
 	/// \note Do not evaluate a 'stub' ever.  It is only used so that combine() can be called
-    c2_binary_function() : c2_function<float_type>(), 
-		Left(*((c2_function<float_type> *)0)), Right(*((c2_function<float_type> *)0)) {}
-	
-	const c2_function<float_type> &Left,  &Right;
-	bool leftown, rightown;
-		
+    c2_binary_function(
+	   float_type (*combiner)(const c2_function<float_type> &left, const c2_function<float_type> &right, 
+							  float_type x, float_type *yprime, float_type *yprime2)
+		) : c2_function<float_type>(), combine(combiner), Left(), Right(), stub(true) { }
+	
+public:
+	float_type (* const combine)(const c2_function<float_type> &left, const c2_function<float_type> &right, 
+								 float_type x, float_type *yprime, float_type *yprime2);
+	
+protected:	
+	const c2_const_ptr<float_type> Left,  Right;
+	/// \brief if true, we don't own any functions, we are just a source of a combining function.
+	bool stub;
+	
 };
 
 /// \brief Create a very lightweight method to return a scalar multiple of another function.
-/// 
-template <typename float_type=double> class c2_scaled_function : public c2_plugin_function<float_type> {
+/// \ingroup containers \ingroup arithmetic_functions \ingroup parametric_functions
+///
+/// The factory function c2_factory::scaled_function() creates *new c2_scaled_function_p
+template <typename float_type=double> class c2_scaled_function_p : public c2_function<float_type> {
 public:
 	/// \brief construct the function with its scale factor.
@@ -456,6 +902,6 @@
 	/// \param outer the function to be scaled
 	/// \param scale the multiplicative scale factor
-	c2_scaled_function(const c2_function<float_type> &outer, float_type scale) : 
-		c2_plugin_function<float_type>(outer), yscale(scale) { }
+	c2_scaled_function_p(const c2_function<float_type> &outer, float_type scale) : 
+		c2_function<float_type>(), func(outer), yscale(scale) { }
 	
 	/// \brief set a new scale factor
@@ -475,20 +921,24 @@
 
 protected:
-    c2_scaled_function<float_type>() {} // hide default constructor, since its use is almost always an error.
+    c2_scaled_function_p<float_type>() : func() {} // hide default constructor, since its use is almost always an error.
+	/// \brief the scaling factor for the function
+	const c2_const_ptr<float_type> func;
 	float_type yscale;
 };
 
 /// \brief A container into which any other c2_function can be dropped.
-///
+/// \ingroup containers
 /// It allows a function to be pre-evaluated at a point, and used at multiple places in an expression
 /// efficiently. If it is re-evaluated at the previous point, it returns the remembered values;
 /// otherwise, it re-evauates the function at the new point.
 ///
-template <typename float_type=double> class c2_cached_function : public c2_plugin_function<float_type> {
+/// The factory function c2_factory::cached_function() creates *new c2_cached_function_p
+template <typename float_type=double> class c2_cached_function_p : public c2_function<float_type> {
 public:
 	/// \brief construct the container
 	///
 	/// \param f the function to be cached
-	c2_cached_function(const c2_function<float_type> &f) : c2_plugin_function<float_type>(f), init(false)  {}
+	c2_cached_function_p(const c2_function<float_type> &f) : c2_function<float_type>(),
+		func(f), init(false)  {}
 	/// \copydoc c2_function::value_with_derivatives
 	///
@@ -508,5 +958,6 @@
 
 protected:
-    c2_cached_function() {} // hide default constructor, since its use is almost always an error.
+    c2_cached_function_p() : func() {} // hide default constructor, since its use is almost always an error.
+	const c2_const_ptr<float_type> func;
 	mutable bool init;
 	mutable float_type x0, y, yp, ypp;
@@ -515,9 +966,9 @@
 
 /// \brief Provides function composition (nesting)
-///
+/// \ingroup arithmetic_functions 
 /// This allows evaluation of \a f(g(x)) where \a f and \a g are c2_function objects.
 ///
-/// \note See c2_binary_function for discussion of ownership.
-template <typename float_type=double> class  c2_composed_function : public c2_binary_function<float_type> {
+/// This should always be constructed using \ref compose_operator "c2_function::operator()"
+template <typename float_type=double> class  c2_composed_function_p : public c2_binary_function<float_type> {
 public:
 	
@@ -526,23 +977,24 @@
 	/// \param outer the outer function 
 	/// \param inner the inner function
-	c2_composed_function(const c2_function<float_type> &outer, const c2_function<float_type> &inner) : c2_binary_function<float_type>(outer, inner) {}
+	c2_composed_function_p(const c2_function<float_type> &outer, const c2_function<float_type> &inner) : 
+		c2_binary_function<float_type>(combine, outer, inner) { this->set_domain(inner.xmin(), inner.xmax()); }
 	/// \brief Create a stub just for the combiner to avoid statics. 
-	c2_composed_function() : c2_binary_function<float_type>() {} ;
-
-	virtual float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, 
-						  float_type x, float_type *yprime, float_type *yprime2) const
+	c2_composed_function_p() : c2_binary_function<float_type>(combine) {} 
+
+	/// \brief execute math necessary to do composition
+	static float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, 
+						  float_type x, float_type *yprime, float_type *yprime2) throw(c2_exception)
     {
-		float_type y0, yp0, ypp0, y1, yp1, ypp1;
-		float_type *yp0p, *ypp0p, *yp1p, *ypp1p;
+		float_type y0, y1;
 		if(yprime || yprime2) {
-			yp0p=&yp0; ypp0p=&ypp0; yp1p=&yp1; ypp1p=&ypp1;
+			float_type yp0, ypp0, yp1, ypp1;
+			y0=right.value_with_derivatives(x, &yp0, &ypp0);
+			y1=left.value_with_derivatives(y0, &yp1, &ypp1);
+			if(yprime) *yprime=yp1*yp0;
+			if(yprime2) *yprime2=ypp0*yp1+yp0*yp0*ypp1;
 		} else {
-			yp0p=ypp0p=yp1p=ypp1p=0;
+			y0=right(x);
+			y1=left(y0);
 		}
-		
-		y0=right.value_with_derivatives(x, yp0p, ypp0p);
-		y1=left.value_with_derivatives(y0, yp1p, ypp1p);
-		if(yprime) *yprime=yp1*yp0;
-		if(yprime2) *yprime2=ypp0*yp1+yp0*yp0*ypp1;
 		return y1;
 	}	
@@ -550,128 +1002,126 @@
 
 /// \brief create a c2_function which is the sum of two other c2_function objects.
-///
-/// \note See c2_binary_function for discussion of ownership.
-template <typename float_type=double> class c2_sum : public c2_binary_function<float_type> {
+/// \ingroup arithmetic_functions
+/// This should always be constructed using c2_function::operator+()
+template <typename float_type=double> class c2_sum_p : public c2_binary_function<float_type> {
 public:	
 	/// \brief construct \a left + \a right
-    /// \note See c2_binary_function for discussion of ownership.
+ 	/// \param left the left function 
+	/// \param right the right function
+	c2_sum_p(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(combine, left, right) {}
+	/// \brief Create a stub just for the combiner to avoid statics. 
+	c2_sum_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics
+
+	/// \brief execute math necessary to do addition
+	static float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, 
+						  float_type x, float_type *yprime, float_type *yprime2) throw(c2_exception)
+    {
+		float_type y0, y1;
+		if(yprime || yprime2) {
+			float_type yp0, ypp0, yp1, ypp1;
+			y0=left.value_with_derivatives(x, &yp0, &ypp0);
+			y1=right.value_with_derivatives(x, &yp1, &ypp1);
+			if(yprime) *yprime=yp0+yp1;
+			if(yprime2) *yprime2=ypp0+ypp1;
+		} else {
+			y0=left(x);
+			y1=right(x);
+		}
+		return y0+y1;
+	}
+};
+
+
+/// \brief create a c2_function which is the difference of two other c2_functions.
+/// \ingroup arithmetic_functions
+/// This should always be constructed using c2_function::operator-()
+template <typename float_type=double> class c2_diff_p : public c2_binary_function<float_type> {
+public:	
+	/// \brief construct \a left - \a right
 	/// \param left the left function 
 	/// \param right the right function
-	c2_sum(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(left, right) {}
+	c2_diff_p(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(combine, left, right) {}
 	/// \brief Create a stub just for the combiner to avoid statics. 
-	c2_sum() : c2_binary_function<float_type>() {} ; // create a stub just for the combiner to avoid statics
-
-	// function to do derivative arithmetic for sums
-	virtual float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, 
-						  float_type x, float_type *yprime, float_type *yprime2) const
+	c2_diff_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics
+
+	/// \brief execute math necessary to do subtraction
+	static float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, 
+							  float_type x, float_type *yprime, float_type *yprime2) throw(c2_exception)
     {
-		float_type y0, yp0, ypp0, y1, yp1, ypp1;
-		float_type *yp0p, *ypp0p, *yp1p, *ypp1p;
+		float_type y0, y1;
 		if(yprime || yprime2) {
-			yp0p=&yp0; ypp0p=& ypp0; yp1p=&yp1; ypp1p=&ypp1;
+			float_type yp0, ypp0, yp1, ypp1;
+			y0=left.value_with_derivatives(x, &yp0, &ypp0);
+			y1=right.value_with_derivatives(x, &yp1, &ypp1);
+			if(yprime) *yprime=yp0-yp1;
+			if(yprime2) *yprime2=ypp0-ypp1;
 		} else {
-			yp0p=ypp0p=yp1p=ypp1p=0;
+			y0=left(x);
+			y1=right(x);
 		}
-		y0=left.value_with_derivatives(x, yp0p, ypp0p);
-		y1=right.value_with_derivatives(x, yp1p, ypp1p);
-		if(yprime) *yprime=yp0+yp1;
-		if(yprime2) *yprime2=ypp0+ypp1;
-		return y0+y1;
+		return y0-y1;
 	}
 };
 
 
-/// \brief create a c2_function which is the difference of two other c2_functions.
-///
-/// \note See c2_binary_function for discussion of ownership.
-template <typename float_type=double> class c2_diff : public c2_binary_function<float_type> {
+/// \brief create a c2_function which is the product of two other c2_functions.
+/// \ingroup arithmetic_functions
+/// This should always be constructed using c2_function::operator*()
+template <typename float_type=double> class c2_product_p : public c2_binary_function<float_type> {
 public:	
-	/// \brief construct \a left - \a right
-    /// \note See c2_binary_function for discussion of ownership.
+	/// \brief construct \a left * \a right
 	/// \param left the left function 
 	/// \param right the right function
-	c2_diff(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(left, right) {}
+	c2_product_p(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(combine, left, right) {}
 	/// \brief Create a stub just for the combiner to avoid statics. 
-	c2_diff() : c2_binary_function<float_type>() {} ; // create a stub just for the combiner to avoid statics
-
-	// function to do derivative arithmetic for diffs
-	virtual float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, 
-							  float_type x, float_type *yprime, float_type *yprime2) const
+	c2_product_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics
+
+	/// \brief execute math necessary to do multiplication
+	static float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, 
+							   float_type x, float_type *yprime, float_type *yprime2)  throw(c2_exception)
     {
-		float_type y0, yp0, ypp0, y1, yp1, ypp1;
-		float_type *yp0p, *ypp0p, *yp1p, *ypp1p;
+		float_type y0, y1;
 		if(yprime || yprime2) {
-			yp0p=&yp0; ypp0p=&ypp0; yp1p=&yp1; ypp1p=&ypp1;
+			float_type yp0, ypp0, yp1, ypp1;
+			y0=left.value_with_derivatives(x, &yp0, &ypp0);
+			y1=right.value_with_derivatives(x, &yp1, &ypp1);
+			if(yprime) *yprime=y1*yp0+y0*yp1;
+			if(yprime2) *yprime2=ypp0*y1+2.0*yp0*yp1+ypp1*y0;
 		} else {
-			yp0p=ypp0p=yp1p=ypp1p=0;
+			y0=left(x);
+			y1=right(x);
 		}
-		y0=left.value_with_derivatives(x, yp0p, ypp0p);
-		y1=right.value_with_derivatives(x, yp1p, ypp1p);
-		if(yprime) *yprime=yp0-yp1;
-		if(yprime2) *yprime2=ypp0-ypp1;
-		return y0-y1;
+		return y0*y1;
 	}
 };
 
 
-/// \brief create a c2_function which is the product of two other c2_functions.
-///
-/// \note See c2_binary_function for discussion of ownership.
-template <typename float_type=double> class c2_product : public c2_binary_function<float_type> {
+/// \brief create a c2_function which is the ratio of two other c2_functions.
+/// \ingroup arithmetic_functions
+/// This should always be constructed using c2_function::operator/()
+template <typename float_type=double> class c2_ratio_p : public c2_binary_function<float_type> {
 public:	
-	/// \brief construct \a left * \a right
-    /// \note See c2_binary_function for discussion of ownership.
+	/// \brief construct \a left / \a right
 	/// \param left the left function 
 	/// \param right the right function
-	c2_product(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(left, right) {}
+	c2_ratio_p(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(combine, left, right) {}
 	/// \brief Create a stub just for the combiner to avoid statics. 
-	c2_product() : c2_binary_function<float_type>() {} ; // create a stub just for the combiner to avoid statics
-
-	virtual float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, 
-							   float_type x, float_type *yprime, float_type *yprime2) const
+	c2_ratio_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics
+	
+	/// \brief execute math necessary to do division
+	static float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, 
+							  float_type x, float_type *yprime, float_type *yprime2) throw(c2_exception)
     {
-		float_type y0, yp0, ypp0, y1, yp1, ypp1;
-		float_type *yp0p, *ypp0p, *yp1p, *ypp1p;
+		float_type y0, y1;
 		if(yprime || yprime2) {
-			yp0p=&yp0; ypp0p=&ypp0; yp1p=&yp1; ypp1p=&ypp1;
+			float_type yp0, ypp0, yp1, ypp1;
+			y0=left.value_with_derivatives(x, &yp0, &ypp0);
+			y1=right.value_with_derivatives(x, &yp1, &ypp1);
+			if(yprime) *yprime=(yp0*y1-y0*yp1)/(y1*y1); // first deriv of ratio
+			if(yprime2) *yprime2=(y1*y1*ypp0+y0*(2*yp1*yp1-y1*ypp1)-2*y1*yp0*yp1)/(y1*y1*y1); 
 		} else {
-			yp0p=ypp0p=yp1p=ypp1p=0;
+			y0=left(x);
+			y1=right(x);
 		}
-		y0=left.value_with_derivatives(x, yp0p, ypp0p);
-		y1=right.value_with_derivatives(x, yp1p, ypp1p);
-		if(yprime) *yprime=y1*yp0+y0*yp1;
-		if(yprime2) *yprime2=ypp0*y1+2.0*yp0*yp1+ypp1*y0;
-		return y0*y1;
-	}
-};
-
-
-/// \brief create a c2_function which is the ratio of two other c2_functions.
-///
-/// \note See c2_binary_function for discussion of ownership.
-template <typename float_type=double> class c2_ratio : public c2_binary_function<float_type> {
-public:	
-	/// \brief construct \a left / \a right
-    /// \note See c2_binary_function for discussion of ownership.
-	/// \param left the left function 
-	/// \param right the right function
-	c2_ratio(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(left, right) {}
-	/// \brief Create a stub just for the combiner to avoid statics. 
-	c2_ratio() : c2_binary_function<float_type>() {} ; // create a stub just for the combiner to avoid statics
-	
-	virtual float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, 
-							  float_type x, float_type *yprime, float_type *yprime2) const
-    {
-		float_type y0, yp0, ypp0, y1, yp1, ypp1;
-		float_type *yp0p, *ypp0p, *yp1p, *ypp1p;
-		if(yprime || yprime2) {
-			yp0p=&yp0; ypp0p=&ypp0; yp1p=&yp1; ypp1p=&ypp1;
-		} else {
-			yp0p=ypp0p=yp1p=ypp1p=0;
-		}
-		y0=left.value_with_derivatives(x, yp0p, ypp0p);
-		y1=right.value_with_derivatives(x, yp1p, ypp1p);
-		if(yprime) *yprime=(yp0*y1-y0*yp1)/(y1*y1); // first deriv of ratio
-		if(yprime2) *yprime2=(y1*y1*ypp0+y0*(2*yp1*yp1-y1*ypp1)-2*y1*yp0*yp1)/(y1*y1*y1); // second deriv of ratio 
 		return y0/y1;
 	}
@@ -679,8 +1129,11 @@
 };
 
-/// \brief a c2_function which is constant : can do interpolating_function  f2=f1 + c2_constant(11.3)
-template <typename float_type> class c2_constant : public c2_function<float_type> {
-public:
-	c2_constant(float_type x=0.0) : c2_function<float_type>(), value(x) {}
+/// \brief a c2_function which is constant 
+/// \ingroup parametric_functions
+///
+/// The factory function c2_factory::constant() creates *new c2_constant_p()
+template <typename float_type> class c2_constant_p : public c2_function<float_type> {
+public:
+	c2_constant_p(float_type x) : c2_function<float_type>(), value(x) {}
 	void reset(float_type val) { value=val; }
 	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) 
@@ -691,7 +1144,206 @@
 };
 
+/// \brief a transformation of a coordinate, including an inverse
+/// \ingroup transforms
+template <typename float_type> class c2_transformation {
+public:
+	/// \brief initialize all our function pointers
+	/// \param transformed true if this function is not the identity
+	/// \param xin input X transform
+	/// \param xinp input X transform derivative
+	/// \param xinpp input X transform second derivative
+	/// \param xout output X transform, which MUST be the inverse of \a xin
+	c2_transformation(bool transformed,
+		   float_type (*xin)(float_type), float_type (*xinp)(float_type), float_type (*xinpp)(float_type), float_type (*xout)(float_type)
+	   ) :
+			fTransformed(transformed), fHasStaticTransforms(true),
+			pIn(xin), pInPrime(xinp), pInDPrime(xinpp), pOut(xout) { }
+
+	/// \brief initialize all our function pointers so that only the (overridden) virtual functions can be called without an error
+	/// \param transformed true if this function is nonlinear
+	c2_transformation(bool transformed) :
+			fTransformed(transformed), fHasStaticTransforms(false),
+			pIn(report_error), pInPrime(report_error), pInDPrime(report_error), pOut(report_error) { }
+	/// \brief the destructor
+	virtual ~c2_transformation() { }
+	/// \brief flag to indicate if this transform is not the identity
+	const bool fTransformed;
+	/// \brief flag to indicate if the static function pointers can be used for efficiency
+	const bool fHasStaticTransforms;
+
+	/// \note the pointers to functions allow highly optimized access when static functions are available.
+	/// They are only used inside value_with_derivatives(), which is assumed to be the most critical routine.
+	/// \brief non-virtual pointer to input X transform
+	float_type (* const pIn)(float_type);
+	/// \brief non-virtual pointer to input X transform derivative
+	float_type (* const pInPrime)(float_type);
+	/// \brief non-virtual pointer to input X transform second derivative
+	float_type (* const pInDPrime)(float_type);
+	/// \brief non-virtual pointer to output X transform 
+	float_type (* const pOut)(float_type);
+	
+	/// \brief virtual input X transform
+	virtual float_type fIn(float_type x) const { return pIn(x); }
+	/// \brief virtual input X transform derivative
+	virtual float_type fInPrime(float_type x) const { return pInPrime(x); }
+	/// \brief virtual input X transform second derivative
+	virtual float_type fInDPrime(float_type x) const { return pInDPrime(x); }
+	/// \brief virtual output X transform 
+	virtual float_type fOut(float_type x) const { return pOut(x); }
+
+protected:
+	/// \brief utility function for unimplemented conversion
+	static float_type report_error(float_type x)  { throw c2_exception("use of improperly constructed axis transform"); return x; }
+	/// \brief utility function f(x)=x useful in axis transforms
+	static float_type ident(float_type x)  { return x; }
+	/// \brief utility function f(x)=1 useful in axis transforms
+	static float_type one(float_type)  { return 1; }
+	/// \brief utility function f(x)=0 useful in axis transforms
+	static float_type zero(float_type)  { return 0; }
+	/// \brief utility function f(x)=1/x useful in axis transforms
+	static float_type recip(float_type x)  { return 1.0/x; }
+	/// \brief utility function f(x)=-1/x**2 useful in axis transforms
+	static float_type recip_prime(float_type x)  { return -1/(x*x); }
+	/// \brief utility function f(x)=2/x**3 useful in axis transforms
+	static float_type recip_prime2(float_type x)  { return 2/(x*x*x); }
+
+};
+
+/// \brief the identity transform
+/// \ingroup transforms
+template <typename float_type> class c2_transformation_linear : public c2_transformation<float_type> {
+public:
+	/// \brief constructor
+	c2_transformation_linear() : c2_transformation<float_type>(false, this->ident, this->one, this->zero, this->ident) { }
+	/// \brief destructor
+	~c2_transformation_linear() { }
+};
+/// \brief log axis transform
+/// \ingroup transforms
+template <typename float_type> class c2_transformation_log : public c2_transformation<float_type> {
+public:
+	/// \brief constructor
+	c2_transformation_log() : c2_transformation<float_type>(true, std::log, this->recip, this->recip_prime, std::exp) { }
+	/// \brief destructor
+	~c2_transformation_log() { }
+};
+/// \brief reciprocal axis transform
+/// \ingroup transforms
+template <typename float_type> class c2_transformation_recip : public c2_transformation<float_type> {
+public:
+	/// \brief constructor
+	c2_transformation_recip() : c2_transformation<float_type>(true, this->recip, this->recip_prime, this->recip_prime2, this->recip) { }
+	/// \brief destructor
+	~c2_transformation_recip() { }
+};
+
+/// \brief a transformation of a function in and out of a coordinate space, using 2 c2_transformations
+///
+/// This class is a container for two axis transforms, but also provides the critical evaluate()
+/// function which converts a result in internal coordinates (with derivatives) into the external representation
+/// \ingroup transforms
+template <typename float_type> 
+	class c2_function_transformation {
+public:
+	/// \brief construct this from two c2_transformation instances
+	/// \param xx the X axis transform
+	/// \param yy the Y axis transform  
+	c2_function_transformation(
+		const c2_transformation<float_type> &xx, const c2_transformation<float_type> &yy) :
+		isIdentity(!(xx.fTransformed || yy.fTransformed)), X(xx), Y(yy) { }
+	/// \brief destructor
+	virtual ~c2_function_transformation() { delete &X; delete &Y; }
+	/// \brief evaluate the transformation from internal coordinates to external coordinates
+	/// \param xraw the value of \a x in external cordinates at which the transform is taking place
+	/// \param y the value of the function in internal coordinates
+	/// \param yp0 the derivative in internal coordinates
+	/// \param ypp0 the second derivative in internal coordinates
+	/// \param [out] yprime pointer to the derivative, or NULL, in external coordinates
+	/// \param [out] yprime2 pointer to the second derivative, or NULL, in external coordinates
+	/// \return the value of the function in external coordinates
+	virtual float_type evaluate(float_type xraw, 
+		float_type y, float_type yp0, float_type ypp0,
+		float_type *yprime, float_type *yprime2) const;
+	/// \brief flag indicating of the transform is the identity, and can be skipped for efficiency
+	const bool isIdentity;
+	/// \brief the X axis transform
+	const c2_transformation<float_type> &X;
+	/// \brief the Y axis transform 
+	const c2_transformation<float_type> &Y;
+};
+
+/// \brief a transformation of a function in and out of lin-lin space
+///
+/// \ingroup transforms
+template <typename float_type> class c2_lin_lin_function_transformation : 
+	public c2_function_transformation<float_type> {
+public:
+	c2_lin_lin_function_transformation() : 
+		c2_function_transformation<float_type>(
+			*new c2_transformation_linear<float_type>, 
+			*new c2_transformation_linear<float_type>
+		) { }
+	virtual ~c2_lin_lin_function_transformation() { }
+};
+
+/// \brief a transformation of a function in and out of log-log space
+///
+/// \ingroup transforms
+template <typename float_type> class c2_log_log_function_transformation : 
+	public c2_function_transformation<float_type> {
+public:
+	c2_log_log_function_transformation() : 
+		c2_function_transformation<float_type>(
+			*new c2_transformation_log<float_type>, 
+			*new c2_transformation_log<float_type>
+		) { }
+	virtual ~c2_log_log_function_transformation() { }
+};
+
+/// \brief a transformation of a function in and out of lin-log space
+///
+/// \ingroup transforms
+template <typename float_type> class c2_lin_log_function_transformation : 
+	public c2_function_transformation<float_type> {
+public:
+	c2_lin_log_function_transformation() : 
+		c2_function_transformation<float_type>(
+			*new c2_transformation_linear<float_type>, 
+			*new c2_transformation_log<float_type>
+		) { }
+	virtual ~c2_lin_log_function_transformation() { }
+};
+
+/// \brief a transformation of a function in and out of log-lin space
+///
+/// \ingroup transforms
+template <typename float_type> class c2_log_lin_function_transformation : 
+	public c2_function_transformation<float_type> {
+public:
+	c2_log_lin_function_transformation() : 
+		c2_function_transformation<float_type>(
+			*new c2_transformation_log<float_type>, 
+			*new c2_transformation_linear<float_type>
+		) { }
+	virtual ~c2_log_lin_function_transformation() { }
+};
+
+/// \brief a transformation of a function in and out of Arrhenuis (1/x vs. log(y)) space
+///
+/// \ingroup transforms
+template <typename float_type> class c2_arrhenius_function_transformation : 
+	public c2_function_transformation<float_type> {
+public:
+	c2_arrhenius_function_transformation() : 
+		c2_function_transformation<float_type>(
+			*new c2_transformation_recip<float_type>, 
+			*new c2_transformation_log<float_type>
+		) { }
+	virtual ~c2_arrhenius_function_transformation() { }
+};
+
 /**
     \brief create a cubic spline interpolation of a set of (x,y) pairs
-
+	\ingroup interpolators
     This is one of the main reasons for c2_function objects to exist.
 
@@ -707,7 +1359,9 @@
     (vanishing second derivative), is created as: \n
     \code 
+	c2_ptr<double> c2p;
+	c2_factory<double> c2;
     std::vector<double> xvals(10), yvals(10); 
     // < fill in xvals and yvals >
-    interpolating_function<double>  myfunc(xvals, yvals);
+    c2p myfunc=c2.interpolating_function().load(xvals, yvals,true,0,true,0);
     // and it can be evaluated at a point for its value only by: 
     double y=myfunc(x); 
@@ -716,11 +1370,23 @@
     double y=myfunc(x,&yprime, &yprime2);
     \endcode
+	
+	The factory function c2_factory::interpolating_function() creates *new interpolating_function_p()
 */
 
-template <typename float_type=double> class interpolating_function  : public c2_function<float_type> {
-public:
-    /// \brief create the most general interpolating_function which defaults to linear-linear space
+template <typename float_type=double> class interpolating_function_p  : public c2_function<float_type> {
+public:
+    /// \brief an empty linear-linear cubic-spline interpolating_function_p
     ///
     /// lots to say here, but see Numerical Recipes for a discussion of cubic splines.
+	///
+	interpolating_function_p() : c2_function<float_type>(), 
+		fTransform(*new  c2_lin_lin_function_transformation<float_type>) { } 
+
+    /// \brief an empty cubic-spline interpolating_function_p with a specific transform
+    ///
+	interpolating_function_p(const c2_function_transformation<float_type> &transform) : c2_function<float_type>(), 
+		fTransform(transform) { } 
+
+    /// \brief do the dirty work of constructing the spline from a function. 
     /// \param x the list of abscissas.  Must be either strictly increasing or strictly decreasing.
 	/// Strictly increasing is preferred, as less memory is used since a copy is not required for the sampling grid.
@@ -730,47 +1396,98 @@
     /// \param upperSlopeNatural if true, set y''(last point)=0, otherwise compute it from \a upperSope
     /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false
-    /// \param inputXConversion a function (not a c2_function) which converts \a x into the internal space. \n
-	/// If this is NULL, use linear space and ignore inputXConversionPrime, inputXConversionDPrime 
-    /// \param inputYConversion a function (not a c2_function) which converts \a y into the internal space. \n
-	/// If this is NULL, use linear space and ignore inputYConversionPrime, inputYConversionDPrime, outputYConversion 
-    /// \param outputYConversion a function (not a c2_function) which converts \a y out of the internal space
-    /// \param inputXConversionPrime the derivative of \a inputXConversion
-    /// \param inputYConversionPrime the derivative of \a inputYConversion
-    /// \param inputXConversionDPrime the second derivative of \a inputXConversion
-    /// \param inputYConversionDPrime the second derivative of \a inputYConversion
-    interpolating_function(const std::vector<float_type> &x, const std::vector<float_type> &f, 
-						  bool lowerSlopeNatural=true, float_type lowerSlope=0.0, 
-						  bool upperSlopeNatural=true, float_type upperSlope=0.0,
-						  float_type (*inputXConversion)(float_type)=0, 
-						  float_type (*inputYConversion)(float_type)=0, 
-						  float_type (*outputYConversion)(float_type)=0, 
-						  float_type (*inputXConversionPrime)(float_type)=0, 
-						  float_type (*inputYConversionPrime)(float_type)=0, 
-						  float_type (*inputXConversionDPrime)(float_type)=0, 
-						  float_type (*inputYConversionDPrime)(float_type)=0 
-						   ) throw(c2_exception) : c2_function<float_type>()
-        { init(x, f, lowerSlopeNatural, lowerSlope, upperSlopeNatural, upperSlope, 
-               inputXConversion, inputYConversion, outputYConversion,
-               inputXConversionPrime, inputYConversionPrime, 
-               inputXConversionDPrime, inputYConversionDPrime 
-               );
-        }
-	
-    /// \brief copy constructor
-    /// \param rhs interpolating_function  to copy from
-	interpolating_function(const interpolating_function <float_type> &rhs) : c2_function<float_type>(rhs), 
-		Xraw(rhs.Xraw), X(rhs.X), F(rhs.F), y2(rhs.y2),
-		fXin(rhs.fXin), fYin(rhs.fYin), fYout(rhs.fYout), 
-		fXinPrime(rhs.fXinPrime), fYinPrime(rhs.fYinPrime), 
-		fXinDPrime(rhs.fXinDPrime), fYinDPrime(rhs.fYinDPrime) ,
-		xInverted(rhs.xInverted), lastKLow(-1)
-	{	set_sampling_grid_pointer(Xraw); }
-	
-	virtual ~interpolating_function() { } // just to suppress warnings about no virtual destructor
+	/// \param splined if true (default), use cubic spline, if false, use linear interpolation.
+	/// \return the same interpolating function, filled
+	interpolating_function_p<float_type> & load(const std::vector<float_type> &x, const std::vector<float_type> &f, 
+			  bool lowerSlopeNatural, float_type lowerSlope, 
+			  bool upperSlopeNatural, float_type upperSlope, bool splined=true
+	) throw(c2_exception);
+
+    /// \brief do the dirty work of constructing the spline from a function. 
+    /// \param data std::vector of std::pairs of x,y.  Will be sorted into x increasing order in place.
+    /// \param lowerSlopeNatural if true, set y''(first point)=0, otherwise compute it from \a lowerSope
+    /// \param lowerSlope derivative of the function at the lower bound, used only if \a lowerSlopeNatural is false
+    /// \param upperSlopeNatural if true, set y''(last point)=0, otherwise compute it from \a upperSope
+    /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false
+	/// \param splined if true (default), use cubic spline, if false, use linear interpolation.
+	/// \return the same interpolating function, filled
+	interpolating_function_p<float_type> & load_pairs(
+			  std::vector<std::pair<float_type, float_type> > &data,  
+			  bool lowerSlopeNatural, float_type lowerSlope, 
+			  bool upperSlopeNatural, float_type upperSlope, bool splined=true
+	) throw(c2_exception);
+
+    /// \brief do the dirty work of constructing the spline from a function.
+    /// \param func a function without any requirement of valid derivatives to sample into an interpolating function.
+	/// Very probably a c2_classic_function.
+    /// \param xmin the lower bound of the region to sample
+    /// \param xmax the upper bound of the region to sample
+	/// \param abs_tol the maximum absolute error permitted when linearly interpolating the points.
+	/// the real error will be much smaller, since this uses cubic splines at the end.
+	/// \param rel_tol the maximum relative error permitted when linearly interpolating the points.
+	/// the real error will be much smaller, since this uses cubic splines at the end.
+    /// \param lowerSlopeNatural if true, set y'(first point) from 3-point parabola, otherwise compute it from \a lowerSope
+    /// \param lowerSlope derivative of the function at the lower bound, used only if \a lowerSlopeNatural is false
+    /// \param upperSlopeNatural if true, set y'(last point) from 3-point parabola, otherwise compute it from \a upperSope
+    /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false
+	/// \return the same interpolating function, filled
+	/// \note If the interpolator being filled has a log vertical axis, put the desired relative error in 
+	/// \a abs_tol, and 0 in \a rel_tol since the absolute error on the log of a function is the relative error
+	/// on the function itself. 
+	interpolating_function_p<float_type> &  sample_function(const c2_function<float_type> &func, 
+			  float_type xmin, float_type xmax, float_type abs_tol, float_type rel_tol,
+			  bool lowerSlopeNatural, float_type lowerSlope, 
+			  bool upperSlopeNatural, float_type upperSlope
+			  ) throw(c2_exception);
+	
+
+	/// \brief initialize from a grid of points and a c2_function (un-normalized) to an 
+	/// interpolator which, when evaluated with a uniform random variate on [0,1] returns random numbers
+	/// distributed as the input function.
+	/// \see  \ref random_subsec "Arbitrary random generation"
+	/// inverse_integrated_density starts with a probability density  std::vector, generates the integral, 
+	/// and generates an interpolating_function_p  of the inverse function which, when evaluated using a uniform random on [0,1] returns values
+	/// with a density distribution equal to the input distribution
+	/// If the data are passed in reverse order (large X first), the integral is carried out from the big end.
+	/// \param bincenters the positions at which to sample the function \a binheights
+	/// \param binheights a function which describes the density of the random number distribution to be produced.
+	/// \return an initialized interpolator, which 
+	/// if evaluated randomly with a uniform variate on [0,1] produces numbers
+	/// distributed according to \a binheights
+	interpolating_function_p<float_type> & load_random_generator_function(
+		const std::vector<float_type> &bincenters, const c2_function<float_type> &binheights)
+		throw(c2_exception);
+
+	/// \brief initialize from a grid of points and an std::vector of probability densities (un-normalized) to an
+	/// interpolator which, when evaluated with a uniform random variate on [0,1] returns random numbers
+	/// distributed as the input histogram.
+	/// \see  \ref random_subsec "Arbitrary random generation"
+	/// inverse_integrated_density starts with a probability density  std::vector, generates the integral, 
+	/// and generates an interpolating_function_p  of the inverse function which, when evaluated using a uniform random on [0,1] returns values
+	/// with a density distribution equal to the input distribution
+	/// If the data are passed in reverse order (large X first), the integral is carried out from the big end.
+	/// \param bins if \a bins .size()==\a binheights .size(), the centers of the bins.  \n
+	/// if \a bins .size()==\a binheights .size()+1, the edges of the bins
+	/// \param binheights a vector which describes the density of the random number distribution to be produced.
+	/// Note density... the numbers in the bins are not counts, but counts/unit bin width.
+	/// \return an initialized interpolator, which 
+	/// if evaluated randomly with a uniform variate on [0,1] produces numbers
+	/// distributed according to \a binheights
+	interpolating_function_p<float_type> & load_random_generator_bins(
+		const std::vector<float_type> &bins, const std::vector<float_type> &binheights)
+		throw(c2_exception);
+	
+	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception);
+	
+	/// \brief destructor
+	virtual ~interpolating_function_p() { delete &fTransform; } 
+	
+	/// \brief create a new, empty interpolating function of this type (virtual constructor)
+	virtual interpolating_function_p<float_type> &clone() const throw(c2_exception)
+	{	return *new interpolating_function_p<float_type>(); }
 	
     /// \brief retrieve copies of the x & y tables from which this was built
     ///
 	/// This is often useful in the creation of new interpolating functions with transformed data.
-	/// The vectorswill have their sizes set correctly on return.
+	/// The vectors will have their sizes set correctly on return.
 	/// \param [in, out] xvals the abscissas 
 	/// \param [in, out] yvals the ordinates
@@ -793,13 +1510,13 @@
 	// preserving the X bounds and mapping functions of the host (left hand) function.
 	
-    /// \brief create a new interpolating_function  which is the \a source 
+    /// \brief create a new interpolating_function_p  which is the \a source 
     /// function applied to every point in the interpolating tables
     ///
     /// This carefully manages the derivative of the composed function at the two ends.
     /// \param source the function to apply
-    /// \return a new interpolating_function  with the same mappings for x and y
-	interpolating_function <float_type> & unary_operator(const c2_function<float_type> &source) const;
-
-    /// \brief create a new interpolating_function  which is the parent interpolating_function  
+    /// \return a new interpolating_function_p  with the same mappings for x and y
+	interpolating_function_p <float_type> & unary_operator(const c2_function<float_type> &source) const;
+
+    /// \brief create a new interpolating_function_p  which is the parent interpolating_function_p  
     /// combined with \a rhs using \a combiner at every point in the interpolating tables
     ///
@@ -807,248 +1524,131 @@
     /// \param rhs the function to apply
     /// \param combining_stub a function which defines which binary operation to use.
-    /// \return a new interpolating_function  with the same mappings for x and y
-	interpolating_function <float_type> & binary_operator(const c2_function<float_type> &rhs,
-           c2_binary_function<float_type> *combining_stub
+    /// \return a new interpolating_function_p  with the same mappings for x and y
+	interpolating_function_p <float_type> & binary_operator(const c2_function<float_type> &rhs,
+           const c2_binary_function<float_type> *combining_stub
            ) const;
-	
-	// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table
-	// when they are applied.  If this is not desired, these operators are not virtual, so the interpolating_function 
-	// can be upcast back to a c2_function to produce unprocessed binaries.
-
-	/// \brief produce a newly resampled interpolating_function  which is the specified sum.
+	/// \brief produce a newly resampled interpolating_function_p  which is the specified sum.
     /// \param rhs the function to add, pointwise
-    /// \return a new interpolating_function 
-    /// \note
-	/// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table
-	/// when they are applied.  If this is not desired, these operators are not virtual, so the interpolating_function 
-	/// can be upcast back to a c2_function to produce unprocessed binaries.
-    interpolating_function <float_type> & operator + (const c2_function<float_type> &rhs) const { 
-		return binary_operator(rhs, new c2_sum<float_type>()); }
-	/// \brief produce a newly resampled interpolating_function  which is the specified difference.
+    /// \return a new interpolating_function_p 
+    interpolating_function_p <float_type> & add_pointwise (const c2_function<float_type> &rhs) const { 
+		return binary_operator(rhs, new c2_sum_p<float_type>()); }
+	/// \brief produce a newly resampled interpolating_function_p  which is the specified difference.
     /// \param rhs the function to subtract, pointwise
-    /// \return a new interpolating_function 
-    /// \note
-	/// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table
-	/// when they are applied.  If this is not desired, these operators are not virtual, so the interpolating_function 
-	/// can be upcast back to a c2_function to produce unprocessed binaries.
-	interpolating_function <float_type> & operator - (const c2_function<float_type> &rhs) const {
-		return binary_operator(rhs, new c2_diff<float_type>()); }
-	/// \brief produce a newly resampled interpolating_function  which is the specified product.
+    /// \return a new interpolating_function_p 
+	interpolating_function_p <float_type> & subtract_pointwise (const c2_function<float_type> &rhs) const {
+		return binary_operator(rhs, new c2_diff_p<float_type>()); }
+	/// \brief produce a newly resampled interpolating_function_p  which is the specified product.
     /// \param rhs the function to multiply, pointwise
-    /// \return a new interpolating_function 
-    /// \note
-	/// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table
-	/// when they are applied.  If this is not desired, these operators are not virtual, so the interpolating_function 
-	/// can be upcast back to a c2_function to produce unprocessed binaries.
-	interpolating_function <float_type> & operator * (const c2_function<float_type> &rhs) const { 
-		return binary_operator(rhs, new c2_product<float_type>()); }
-	/// \brief produce a newly resampled interpolating_function  which is the specified ratio.
+    /// \return a new interpolating_function_p 
+	interpolating_function_p <float_type> & multiply_pointwise (const c2_function<float_type> &rhs) const { 
+		return binary_operator(rhs, new c2_product_p<float_type>()); }
+	/// \brief produce a newly resampled interpolating_function_p  which is the specified ratio.
     /// \param rhs the function to divide, pointwise
-    /// \return a new interpolating_function 
-    /// \note
-	/// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table
-	/// when they are applied.  If this is not desired, these operators are not virtual, so the interpolating_function 
-	/// can be upcast back to a c2_function to produce unprocessed binaries.
-	interpolating_function <float_type> & operator / (const c2_function<float_type> &rhs) const { 
-		return binary_operator(rhs, new c2_ratio<float_type>()); }
-	/// \brief produce a newly resampled interpolating_function  which is the specified sum.
-    /// \param rhs a constant to add, pointwise
-    /// \return a new interpolating_function 
-    /// \note
-	/// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table
-	/// when they are applied.  If this is not desired, these operators are not virtual, so the interpolating_function 
-	/// can be upcast back to a c2_function to produce unprocessed binaries.
-	interpolating_function <float_type> & operator + (float_type rhs) const { return (*this)+c2_constant<float_type>(rhs); }
-	/// \brief produce a newly resampled interpolating_function  which is the specified difference.
-    /// \param rhs a constant to subtract, pointwise
-    /// \return a new interpolating_function 
-    /// \note
-	/// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table
-	/// when they are applied.  If this is not desired, these operators are not virtual, so the interpolating_function 
-	/// can be upcast back to a c2_function to produce unprocessed binaries.
-	interpolating_function <float_type> & operator - (float_type rhs) const { return (*this)-c2_constant<float_type>(rhs); }
-	/// \brief produce a newly resampled interpolating_function  which is the specified product.
-    /// \param rhs a constant to multiply, pointwise
-    /// \return a new interpolating_function 
-    /// \note
-	/// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table
-	/// when they are applied.  If this is not desired, these operators are not virtual, so the interpolating_function 
-	/// can be upcast back to a c2_function to produce unprocessed binaries.
-	interpolating_function <float_type> & operator * (float_type rhs) const { return (*this)*c2_constant<float_type>(rhs); }
-	/// \brief produce a newly resampled interpolating_function  which is the specified ratio.
-    /// \param rhs a constant to divide, pointwise
-    /// \return a new interpolating_function 
-    /// \note
-	/// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table
-	/// when they are applied.  If this is not desired, these operators are not virtual, so the interpolating_function 
-	/// can be upcast back to a c2_function to produce unprocessed binaries.
-	interpolating_function <float_type> & operator / (float_type rhs) const { return (*this)/c2_constant<float_type>(rhs); }
-        
-	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception);
-	
-	/// \brief move value & derivatives into our internal coordinates (use splint to go the other way!)
-    /// \note why?
-	void localize_derivatives(float_type xraw, float_type y, float_type yprime, float_type yprime2, float_type *y0, float_type *yp0, float_type *ypp0) const;
-	
-protected:
-		
-	interpolating_function() : c2_function<float_type>() { } // default constructor is never used, prevent accidents by protecting it.
-    
-    /// \brief do the dirty work of constructing the spline.  See interpolating_function  constructor for details.
-	void init(const std::vector<float_type> &, const std::vector<float_type> &, 
+    /// \return a new interpolating_function_p 
+	interpolating_function_p <float_type> & divide_pointwise (const c2_function<float_type> &rhs) const { 
+		return binary_operator(rhs, new c2_ratio_p<float_type>()); }
+    /// \brief copy data from another interpolating function.  This only makes sense if the source
+	/// function has the same transforms as the destination.
+    /// \param rhs interpolating_function_p  to copy from
+	void clone_data(const interpolating_function_p <float_type> &rhs) {
+		Xraw=rhs.Xraw; X=rhs.X; F=rhs.F; y2=rhs.y2;
+		set_sampling_grid_pointer(Xraw);
+	}
+
+	const c2_function_transformation<float_type> &fTransform;
+	
+protected:    
+    /// \brief create the spline coefficients
+	void spline(
 			  bool lowerSlopeNatural, float_type lowerSlope, 
-			  bool upperSlopeNatural, float_type upperSlope,
-			  float_type (*inputXConversion)(float_type)=0, 
-			  float_type (*inputXConversionPrime)(float_type)=0, 
-			  float_type (*inputXConversionDPrime)(float_type)=0, 
-			  float_type (*inputYConversion)(float_type)=0, 
-			  float_type (*inputYConversionPrime)(float_type)=0, 
-			  float_type (*inputYConversionDPrime)(float_type)=0, 
-			  float_type (*outputYConversion)(float_type)=0
-			  ) throw(c2_exception) ;
+			  bool upperSlopeNatural, float_type upperSlope
+			  ) throw(c2_exception);
+
+	// This is for sorting the data. It must be static if it's going to be a class member.
+	static bool comp_pair(std::pair<float_type,float_type> const &i, std::pair<float_type,float_type> const &j) {return i.first<j.first;}
 
     std::vector<float_type> Xraw, X, F, y2;
-	
-	float_type (*fXin)(float_type), (*fYin)(float_type), (*fYout)(float_type);
-	float_type (*fXinPrime)(float_type), (*fYinPrime)(float_type);
-	float_type (*fXinDPrime)(float_type), (*fYinDPrime)(float_type);
-	
-	int xInverted;
-    mutable int lastKLow;
-};
-
-/// \brief An interpolatingFunction with X transformed into log space.  
-///
+	c2_const_ptr<float_type> sampler_function;
+	bool xInverted;
+    mutable size_t lastKLow;
+};
+
+/// \brief A spline with X transformed into log space.  
+/// \ingroup interpolators
 /// Most useful for functions looking like y=log(x) or any other function with a huge X dynamic range,
 /// and a slowly varying Y.
-template <typename float_type=double> class log_lin_interpolating_function : public interpolating_function <float_type> {
-public:
-    /// \brief Construct the function.
-    /// \param x the list of abscissas.  Must be either strictly increasing or strictly decreasing.
-	/// Strictly increasing is preferred, as less memory is used since a copy is not required for the sampling grid.
-    /// \param f the list of function values.
-    /// \param lowerSlopeNatural if true, set y''(first point)=0 in LogLin space, otherwise compute it from \a lowerSope
-    /// \param lowerSlope derivative of the function at the lower bound, used only if \a lowerSlopeNatural is false
-    /// \param upperSlopeNatural if true, set y''(last point)=0 in LogLin space, otherwise compute it from \a upperSope
-    /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false
-    log_lin_interpolating_function(const std::vector<float_type> &x, const std::vector<float_type> &f, 
-								bool lowerSlopeNatural=true, float_type lowerSlope=0.0, 
-								bool upperSlopeNatural=true, float_type upperSlope=0.0);
-protected:
-	log_lin_interpolating_function() {} // do not allow naked construction... it is usually an accident.
-};
-
-
-/// \brief An interpolatingFunction with Y transformed into log space.  
-///
+///
+///	The factory function c2_factory::log_lin_interpolating_function() creates *new log_lin_interpolating_function_p()
+template <typename float_type=double> class log_lin_interpolating_function_p : public interpolating_function_p <float_type> {
+public:
+    /// \brief an empty log-linear cubic-spline interpolating_function_p
+    ///
+	log_lin_interpolating_function_p() : interpolating_function_p<float_type>(*new c2_log_lin_function_transformation<float_type>)
+		{ }
+	/// \brief create a new, empty interpolating function of this type (virtual constructor)
+	virtual interpolating_function_p<float_type> &clone() const throw(c2_exception)
+		{	return *new log_lin_interpolating_function_p<float_type>(); }
+};
+
+
+/// \brief A spline with Y transformed into log space.  
+/// \ingroup interpolators
 /// Most useful for functions looking like y=exp(x)
-template <typename float_type=double> class lin_log_interpolating_function : public interpolating_function <float_type> {
-public:
-    /// \brief Construct the function.
-    /// \param x the list of abscissas.  Must be either strictly increasing or strictly decreasing.
-	/// Strictly increasing is preferred, as less memory is used since a copy is not required for the sampling grid.
-    /// \param f the list of function values.
-    /// \param lowerSlopeNatural if true, set y''(first point)=0 in LinLog space, otherwise compute it from \a lowerSope
-    /// \param lowerSlope derivative of the function at the lower bound, used only if \a lowerSlopeNatural is false
-    /// \param upperSlopeNatural if true, set y''(last point)=0 in LinLog space, otherwise compute it from \a upperSope
-    /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false
-    lin_log_interpolating_function(const std::vector<float_type> &x, const std::vector<float_type> &f, 
-								bool lowerSlopeNatural=true, float_type lowerSlope=0.0, 
-								bool upperSlopeNatural=true, float_type upperSlope=0.0);
-protected:
-	lin_log_interpolating_function() {} // do not allow naked construction... it is usually an accident.
-};
-
-
-/// \brief An interpolatingFunction with X and Y transformed into log space.  
-///
+///
+///	The factory function c2_factory::lin_log_interpolating_function() creates *new lin_log_interpolating_function_p()
+template <typename float_type=double> class lin_log_interpolating_function_p : public interpolating_function_p <float_type> {
+public:
+    /// \brief an empty linear-log cubic-spline interpolating_function_p
+    ///
+	lin_log_interpolating_function_p() : interpolating_function_p<float_type>(*new c2_lin_log_function_transformation<float_type>)
+		{ } 
+	/// \brief create a new, empty interpolating function of this type (virtual constructor)
+	virtual interpolating_function_p<float_type> &clone() const throw(c2_exception)
+		{	return *new lin_log_interpolating_function_p<float_type>(); }
+};
+
+
+/// \brief A spline with X and Y transformed into log space.  
+/// \ingroup interpolators
 /// Most useful for functions looking like y=x^n or any other function with a huge X and Y dynamic range.
-template <typename float_type=double> class log_log_interpolating_function : public interpolating_function <float_type> {
-public:
-    /// \brief Construct the function.
-    /// \param x the list of abscissas.  Must be either strictly increasing or strictly decreasing.
-	/// Strictly increasing is preferred, as less memory is used since a copy is not required for the sampling grid.
-    /// \param f the list of function values.
-    /// \param lowerSlopeNatural if true, set y''(first point)=0 in LogLog space, otherwise compute it from \a lowerSope
-    /// \param lowerSlope derivative of the function at the lower bound, used only if \a lowerSlopeNatural is false
-    /// \param upperSlopeNatural if true, set y''(last point)=0 in LogLog space, otherwise compute it from \a upperSope
-    /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false
-    log_log_interpolating_function(const std::vector<float_type> &x, const std::vector<float_type> &f, 
-								bool lowerSlopeNatural=true, float_type lowerSlope=0.0, 
-								bool upperSlopeNatural=true, float_type upperSlope=0.0);
-protected:
-	log_log_interpolating_function() {} // do not allow naked construction... it is usually an accident.
-};
-
-
-/// \brief An interpolating_function  with X in reciprocal space and Y transformed in log space.  
-///
+///
+///	The factory function c2_factory::log_log_interpolating_function() creates *new log_log_interpolating_function_p()
+template <typename float_type=double> class log_log_interpolating_function_p : public interpolating_function_p <float_type> {
+public:	
+    /// \brief an empty log-log cubic-spline interpolating_function_p
+    ///
+	log_log_interpolating_function_p() : interpolating_function_p<float_type>(*new c2_log_log_function_transformation<float_type>)
+		{ } 
+	/// \brief create a new, empty interpolating function of this type (virtual constructor)
+	virtual interpolating_function_p<float_type> &clone() const throw(c2_exception)
+	{	return *new log_log_interpolating_function_p<float_type>(); }
+};
+
+
+/// \brief A spline with X in reciprocal space and Y transformed in log space.  
+/// \ingroup interpolators
 /// Most useful for thermodynamic types of data where Y is roughly A*exp(-B/x). 
 /// Typical examples are reaction rate data, and thermistor calibration data.
-template <typename float_type=double> class arrhenius_interpolating_function : public interpolating_function <float_type> {
-public:
-    /// \brief Construct the function.
-    /// \param x the list of abscissas.  Must be either strictly increasing or strictly decreasing.
-	/// Strictly increasing is preferred, as less memory is used since a copy is not required for the sampling grid.
-    /// \param f the list of function values.
-    /// \param lowerSlopeNatural if true, set y''(first point)=0 in Arrhenius space, otherwise compute it from \a lowerSope
-    /// \param lowerSlope derivative of the function at the lower bound, used only if \a lowerSlopeNatural is false
-    /// \param upperSlopeNatural if true, set y''(last point)=0 in Arrhenius space, otherwise compute it from \a upperSope
-    /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false
-    arrhenius_interpolating_function(const std::vector<float_type> &x, const std::vector<float_type> &f, 
-								bool lowerSlopeNatural=true, float_type lowerSlope=0.0, 
-								bool upperSlopeNatural=true, float_type upperSlope=0.0);
-protected:
-	arrhenius_interpolating_function() {} // do not allow naked construction... it is usually an accident.
-};
-
-/**
- \brief create a linear-linear interpolating grid with both x & y set to 
- (xmin, xmin+dx, ... xmin + (count-1)*dx )
-
- very useful for transformaiton with other functions e.g. 
- \code
- f=c2_sin<double>::sin(LinearInterpolatingGrid(-0.1,0.1, 65)) 
- \endcode
- creates a spline table of sin(x) slightly beyond the first period
- \param xmin the starting point for the grid
- \param dx the step size for the grid
- \param count the number of points in the  grid
- \return an identity interpolating_function  with the requested grid 
- */
-template <typename float_type> interpolating_function <float_type> &linear_interpolating_grid(float_type xmin, float_type dx, int count) {
-	std::vector<float_type> x(count);
-	for(int i=0; i<count; i++) x[i]=xmin + i * dx;
-	return *new interpolating_function <float_type>(x,x);
-}
-
-/**
- \brief create a log-log interpolating grid with both x & y set to 
- (xmin, xmin*dx, ... xmin * dx^(count-1) )
-
- very useful for transformaiton with other functions e.g. 
- \code
- f=c2_log<double>::log(LogLogInterpolatingGrid(2, 1.1, 65)) 
- \endcode
- creates a spline table of log(x)
- \param xmin the starting point for the grid
- \param dx the ratio between points
- \param count the number of points in the  grid
- \return an identity log_log_interpolating_function with the requested grid 
-*/
-template <typename float_type> log_log_interpolating_function <float_type> &log_log_interpolating_grid(float_type xmin, float_type dx, int count) {
-	std::vector<float_type> x(count);
-	x[0]=xmin;
-	for(int i=1; i<count; i++) x[i]=dx*x[i-1];
-	return *new log_log_interpolating_function<float_type>(x,x);
-}
+///
+///	The factory function c2_factory::arrhenius_interpolating_function() creates *new arrhenius_interpolating_function_p()
+template <typename float_type=double> class arrhenius_interpolating_function_p : public interpolating_function_p <float_type> {
+public:
+    /// \brief an empty arrhenius cubic-spline interpolating_function_p
+    ///
+	arrhenius_interpolating_function_p() : interpolating_function_p<float_type>(*new c2_arrhenius_function_transformation<float_type>)
+		{ } 
+	/// \brief create a new, empty interpolating function of this type (virtual constructor)
+	virtual interpolating_function_p<float_type> &clone() const throw(c2_exception)
+	{	return *new arrhenius_interpolating_function_p<float_type>(); }	
+};
 
 /// \brief compute sin(x) with its derivatives.
-///
-/// Creates a singleton instance c2_sin::sin of itself for convenient access.
-template <typename float_type=double> class c2_sin : public c2_function<float_type> {
-public:
-	/// \brief constructor.  There is alread a singleton c2_sin::sin, which usually suffices.
-	c2_sin() {}
+/// \ingroup math_functions
+///
+/// The factory function c2_factory::sin() creates *new c2_sin_p
+template <typename float_type=double> class c2_sin_p : public c2_function<float_type> {
+public:
+	/// \brief constructor.
+	c2_sin_p() : c2_function<float_type>() {}
 	
 	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
@@ -1058,29 +1658,29 @@
     /// \param xmin the lower bound for the grid
     /// \param xmax upper bound for the grid
-    /// \return a new sampling grid.
-	virtual std::vector<float_type> &get_sampling_grid(float_type xmin, float_type xmax); 
-    /// \brief the static singleton 
-	static const c2_sin sin;
-};
+    /// \param [in, out] grid the sampling grid.
+	virtual void get_sampling_grid(float_type xmin, float_type xmax,  std::vector<float_type> &grid) const; 
+};
+
 /// \brief compute cos(x) with its derivatives.
-///
-/// Creates a singleton instance c2_cos::cos of itself for convenient access.
-template <typename float_type=double> class c2_cos : public c2_sin<float_type> {
-public:
-	/// \brief constructor.  There is already a singleton c2_cos::cos, which usually suffices.
-	c2_cos() {}
+/// \ingroup math_functions
+///
+/// The factory function c2_factory::cos() creates *new c2_cos_p
+template <typename float_type=double> class c2_cos_p : public c2_sin_p<float_type> {
+public:
+	/// \brief constructor.
+	c2_cos_p() : c2_sin_p<float_type>() {}
 	
 	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
 	{ float_type q=std::cos(x); if(yprime) *yprime=-std::sin(x); if(yprime2) *yprime2=-q; return q; }	
-    /// \brief the static singleton 
-	static const c2_cos cos;
-};
+};
+
 /// \brief compute tan(x) with its derivatives.
-///
-/// Creates a singleton instance c2_tan::tan of itself for convenient access.
-template <typename float_type=double> class c2_tan : public c2_function<float_type> {
-public:
-	/// \brief constructor.  There is already a singleton c2_tan::tan, which usually suffices.
-	c2_tan() {}
+/// \ingroup math_functions
+///
+/// The factory function c2_factory::tan() creates *new c2_tan_p
+template <typename float_type=double> class c2_tan_p : public c2_function<float_type> {
+public:
+	/// \brief constructor.
+	c2_tan_p() : c2_function<float_type>() {}
 	
 	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
@@ -1092,53 +1692,53 @@
 		return t; 
 	}	
-    /// \brief the static singleton 
-	static const c2_tan tan;
-};
+};
+
 /// \brief compute log(x) with its derivatives.
-///
-/// Creates a singleton instance c2_log::log of itself for convenient access.
-template <typename float_type=double> class c2_log : public c2_function<float_type> {
-public:
-	/// \brief constructor.  There is already a singleton c2_log::log, which usually suffices.
-	c2_log() {}
+/// \ingroup math_functions
+///
+/// The factory function c2_factory::log() creates *new c2_log_p
+template <typename float_type=double> class c2_log_p : public c2_function<float_type> {
+public:
+	/// \brief constructor.
+	c2_log_p() : c2_function<float_type>() {}
 
 	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
 	{ if(yprime) *yprime=1.0/x; if(yprime2) *yprime2=-1.0/(x*x); return std::log(x); }	
-    /// \brief the static singleton 
-	static const c2_log log;
-};
+};
+
 /// \brief compute exp(x) with its derivatives.
-///
-/// Creates a singleton instance c2_exp::exp of itself for convenient access.
-template <typename float_type=double>  class c2_exp : public c2_function<float_type> {
-public:
-	/// \brief constructor.  There is already a singleton c2_exp::exp, which usually suffices.
-	c2_exp() {}
+/// \ingroup math_functions
+///
+/// The factory function c2_factory::exp() creates *new c2_exp_p
+template <typename float_type=double>  class c2_exp_p : public c2_function<float_type> {
+public:
+	/// \brief constructor.
+	c2_exp_p() : c2_function<float_type>() {}
 
 	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
 	{ float_type q=std::exp(x); if(yprime) *yprime=q; if(yprime2) *yprime2=q; return q; }
-    /// \ brief the static singleton 
-	static const c2_exp exp;
-};
+};
+
 /// \brief compute sqrt(x) with its derivatives.
-///
-/// Creates a singleton instance c2_sqrt::sqrt of itself for convenient access.
-template <typename float_type=double> class c2_sqrt : public c2_function<float_type> {
-public:
-	/// \brief constructor.  There is already a singleton c2_sqrt::sqrt, which usually suffices.
-	c2_sqrt() {}
+/// \ingroup math_functions
+///
+/// The factory function c2_factory::sqrt() creates *new c2_sqrt_p()
+template <typename float_type=double> class c2_sqrt_p : public c2_function<float_type> {
+public:
+	/// \brief constructor.
+	c2_sqrt_p() : c2_function<float_type>() {}
 	
 	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
 	{ float_type q=std::sqrt(x); if(yprime) *yprime=0.5/q; if(yprime2) *yprime2=-0.25/(x*q); return q; }
-    /// \brief the static singleton 
-	static const c2_sqrt sqrt;
-};
+};
+
 /// \brief compute scale/x with its derivatives.
-///
-/// Creates a singleton instance c2_recip:recip of itself for convenient access.
-template <typename float_type=double> class c2_recip : public c2_function<float_type> {
-public:
-	/// \brief constructor.  There is already a singleton c2_recip::recip, which usually suffices.
-	c2_recip(float_type scale) : rscale(scale) {}
+/// \ingroup parametric_functions
+///
+/// The factory function c2_factory::recip() creates *new c2_recip_p
+template <typename float_type=double> class c2_recip_p : public c2_function<float_type> {
+public:
+	/// \brief constructor.
+	c2_recip_p(float_type scale) : c2_function<float_type>(), rscale(scale) {}
 	
 	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
@@ -1153,34 +1753,33 @@
 	/// \param scale the new numerator
 	void reset(float_type scale) { rscale=scale; } 
-    /// \brief the static singleton 
-	static const c2_recip recip;
 private:
 	float_type rscale;
 };
+
 /// \brief compute x with its derivatives.
-///
-/// Creates a singleton instance c2_identity::identity of itself for convenient access.
-template <typename float_type=double> class c2_identity : public c2_function<float_type> {
-public:
-	/// \brief constructor.  There is already a singleton c2_identity::identity, which usually suffices.
-	c2_identity() {}
+/// \ingroup math_functions
+///
+/// The factory function c2_factory::identity() creates *new c2_identity_p
+template <typename float_type=double> class c2_identity_p : public c2_function<float_type> {
+public:
+	/// \brief constructor.
+	c2_identity_p() : c2_function<float_type>() {}
 	
 	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
 	{ if(yprime) *yprime=1.0; if(yprime2) *yprime2=0; return x; }
-    /// \brief the static singleton 
-	static const c2_identity identity;
 };
 
 /**
  \brief create a linear mapping of another function 
- 
+ \ingroup parametric_functions
  for example, given a c2_function \a f 
  \code 
- c2_linear<double> L(1.2, 2.0, 3.0);
- c2_composed_function<double> &F=L(f); 
+ c2_function<double> &F=c2_linear<double>(1.2, 2.0, 3.0)(f);
  \endcode
  produces a new c2_function F=2.0+3.0*(\a f - 1.2) 
+ 
+ The factory function c2_factory::linear() creates *new c2_linear_p
 */
-template <typename float_type=double> class c2_linear : public c2_function<float_type> {
+template <typename float_type=double> class c2_linear_p : public c2_function<float_type> {
 public:
     /// \brief Construct the operator f=y0 + slope * (x-x0)
@@ -1188,5 +1787,6 @@
     /// \param y0 the y-intercept i.e. f(x0)
     /// \param slope the slope of the mapping
-	c2_linear(float_type x0, float_type y0, float_type slope) : xint(x0), intercept(y0), m(slope) {}
+	c2_linear_p(float_type x0, float_type y0, float_type slope) : 
+		c2_function<float_type>(), xint(x0), intercept(y0), m(slope) {}
     /// \brief Change the slope and intercepts after construction.
 	/// \param x0 the x offset
@@ -1200,14 +1800,13 @@
 		float_type xint, intercept, m;
 protected:
-		c2_linear() {} // do not allow naked construction... it is usually an accident.
+		c2_linear_p() {} // do not allow naked construction... it is usually an accident.
 };
 
 /**
 \brief create a quadratic mapping of another function 
- 
+ \ingroup parametric_functions
  for example, given a c2_function \a f 
  \code 
- c2_quadratic<double> Q(1.2, 2.0, 3.0, 4.0);
- c2_composed_function<double> &F=Q(f); 
+ c2_function<double> &F=c2_quadratic<double>(1.2, 2.0, 3.0, 4.0)(f);
  \endcode
  produces a new c2_function F=2.0 + 3.0*(f-1.2) + 4.0*(f-1.2)^2 
@@ -1215,6 +1814,7 @@
  note that the parameters are overdetermined, but allows the flexibility of two different representations
 
+ The factory function c2_factory::quadratic() creates *new c2_quadratic_p
  */
-template <typename float_type=double> class c2_quadratic : public c2_function<float_type> {
+template <typename float_type=double> class c2_quadratic_p : public c2_function<float_type> {
 public:
     /// \brief Construct the operator
@@ -1223,6 +1823,7 @@
     /// \param xcoef the scale on the (\a x - \a x0) term
     /// \param x2coef the scale on the (\a x - \a x0)^2 term
-	c2_quadratic(float_type x0, float_type y0, float_type xcoef, float_type x2coef) : intercept(y0), center(x0), a(x2coef), b(xcoef) {}
-    /// Modify the mapping after construction
+	c2_quadratic_p(float_type x0, float_type y0, float_type xcoef, float_type x2coef) : 
+		c2_function<float_type>(), intercept(y0), center(x0), a(x2coef), b(xcoef) {}
+    /// \brief Modify the coefficients after construction
     /// \param x0 the new center around which the powers are computed
     /// \param y0 the new value of the function at \a x = \a x0
@@ -1236,24 +1837,26 @@
 		float_type intercept, center, a, b;
 protected:
-		c2_quadratic() {} // do not allow naked construction... it is usually an accident.
+		c2_quadratic_p() {} // do not allow naked construction... it is usually an accident.
 };
 
 /**
 \brief create a power law mapping of another function 
- 
+ \ingroup parametric_functions
  for example, given a c2_function \a f 
  \code 
- c2_power_law<double> PLaw(1.2, 2.5);
- c2_composed_function<double> &F=PLaw(f); 
+ c2_power_law_p<double> PLaw(1.2, 2.5);
+ c2_composed_function_p<double> &F=PLaw(f); 
  \endcode
  produces a new c2_function F=1.2 * f^2.5 
  
+ The factory function c2_factory::power_law() creates *new c2_power_law_p
  */
-template <typename float_type=double> class c2_power_law : public c2_function<float_type> {
+template <typename float_type=double> class c2_power_law_p : public c2_function<float_type> {
 public:
     /// \brief Construct the operator
     /// \param scale the multipler
     /// \param power the exponent
-	c2_power_law(float_type scale, float_type power) : a(scale), b(power) {}
+	c2_power_law_p(float_type scale, float_type power) :
+		c2_function<float_type>(), a(scale), b(power) {}
     /// \brief Modify the mapping after construction
     /// \param scale the new multipler
@@ -1266,10 +1869,10 @@
 		float_type a, b;
 protected:
-		c2_power_law() {} // do not allow naked construction... it is usually an accident.
+		c2_power_law_p() {} // do not allow naked construction... it is usually an accident.
 };
 
 /**
 \brief create the formal inverse function of another function 
- 
+ \ingroup containers
  for example, given a c2_function \a f 
  \code 
@@ -1282,12 +1885,14 @@
  derivatives, too, unlike the case of a simple root-finding inverse.  This means
  it can be integrated (for example) quite efficiently.
- 
- Note that it is a subclass of c2_scaled_function only to manage ownership of another c2_function.
- */
-template <typename float_type=double> class c2_inverse_function : public c2_plugin_function<float_type> {
+
+ \see \ref combined_inversion_hinting_sampling
+
+ The factory function c2_factory::inverse_function() creates *new c2_inverse_function_p
+*/
+template <typename float_type=double> class c2_inverse_function_p : public c2_function<float_type> {
 public:
     /// \brief Construct the operator
     /// \param source the function to be inverted
-	c2_inverse_function(const c2_function<float_type> &source);  
+	c2_inverse_function_p(const c2_function<float_type> &source);  
     virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception);
     
@@ -1301,15 +1906,43 @@
     ///  It is used in value_with_derivatives() to guess where to start the root finder.
     /// \param x the abscissa for which an estimate is needed
-    virtual float_type get_start_hint(float_type x) const { return start_hint; } 
-    
+    virtual float_type get_start_hint(float_type x) const 
+		{ return hinting_function.valid()? hinting_function(x) : start_hint; } 
+
+	/// \brief set or unset the approximate function used to start the root finder
+	/// \anchor set_hinting_function_discussion
+	/// A hinting function is mostly useful if the evaluation of this inverse is
+	/// going to be carried out in very non-local order, so the root finder has to start over 
+	/// for each step.  If most evaluations are going to be made in fairly localized clusters (scanning
+	/// through the function, for example), the default mechanism used (which just remembers the last point)
+	/// is almost certainly faster.
+	/// 
+	/// Typically, the hinting function is likely to be set up by creating the inverse function,
+	/// and then adaptively sampling an interpolating function from it, and then using the result
+	/// to hint it.  Another way, if the parent function is already an interpolating function, is just to create
+	/// a version of the parent with the x & y coordinates reversed.
+	/// 
+	/// \see \ref combined_inversion_hinting_sampling
+	///
+	/// \param hint_func the function that is an approximate inverse of the parent of this inverse_function
+	void set_hinting_function(const c2_function<float_type> *hint_func) 
+		{ hinting_function.set_function(hint_func); }
+    /// \brief set the hinting function from a pointer.
+	/// 
+	/// See \ref set_hinting_function_discussion "discussion"
+	/// \param hint_func the container holding the function
+	void set_hinting_function(const c2_const_ptr<float_type> hint_func) 
+		{ hinting_function=hint_func; }
+	
 protected:
-    c2_inverse_function() {} // do not allow naked construction... it is usually an accident.
-    mutable float_type start_hint;
+	c2_inverse_function_p() {} // do not allow naked construction... it is usually an accident.
+	mutable float_type start_hint;
+	const c2_const_ptr<float_type> func;
+	c2_const_ptr<float_type> hinting_function;
 };
 
 /** 
   \brief
-    An interpolating_function  which is the cumulative integral of a histogram.
-
+    An interpolating_function_p  which is the cumulative integral of a histogram.
+	\ingroup interpolators
     Note than binedges should be one element longer than binheights, since the lower & upper edges are specified. 
     Note that this is a malformed spline, since the second derivatives are all zero, so it has less continuity.
@@ -1318,5 +1951,5 @@
 */
 
-template <typename float_type=double>  class accumulated_histogram : public interpolating_function <float_type> {
+template <typename float_type=double>  class accumulated_histogram : public interpolating_function_p <float_type> {
 public:
     /// \brief Construct the integrated histogram
@@ -1332,79 +1965,50 @@
 
 /**
-  \brief Construct a function useful for generation of random numbers from the given distribution
-  
-  inverse_integrated_density<InterpolatingFunctionFlavor>() starts with a probability density c2_function, generates the integral, 
-  and generates an interpolating_function  which, when evaluated using a uniform random on [0,1] returns values
-  with a density distribution equal to the input distribution
-  If the data are passed in reverse order (large X first), the integral is carried out from the big end.
+	\anchor inverse_integrated_density_bins
+	\brief construct from a grid of points and an std::vector of probability densities (un-normalized)
+	\see  \ref random_subsec "Arbitrary random generation"
+	\ingroup interpolators
+	 inverse_integrated_density starts with a probability density  std::vector, generates the integral, 
+	 and generates an interpolating_function_p  of the inverse function which, when evaluated using a uniform random on [0,1] returns values
+	 with a density distribution equal to the input distribution
+	 If the data are passed in reverse order (large X first), the integral is carried out from the big end.
+	 
+	\param bins if \a bins .size()==\a binheights .size(), the centers of the bins.  \n
+	if \a bins .size()==\a binheights .size()+1, the edges of the bins
+	\param binheights a vector which describes the density of the random number distribution to be produced.
+	Note density... the numbers in the bins are not counts, but counts/unit bin width.
+	\return an interpolating_function_p  of the type requested in the template which,
+	if evaluated randomly with a uniform variate on [0,1] produces numbers
+	distributed according to \a binheights
+*/
+
+template <typename float_type, typename Final> 
+interpolating_function_p<float_type> & inverse_integrated_density_bins(
+		const std::vector<float_type> &bins, const std::vector<float_type> &binheights) 
+		throw(c2_exception);
+
+/**
+\anchor inverse_integrated_density_function
+	\brief construct from a grid of points and a c2_function of probability densities (un-normalized)
+	\see  \ref random_subsec "Arbitrary random generation"
+	\ingroup interpolators 
+ inverse_integrated_density starts with a probability density  std::vector, generates the integral, 
+ and generates an interpolating_function_p  of the inverse function which, when evaluated using a uniform random on [0,1] returns values
+ with a density distribution equal to the input distribution
+ If the data are passed in reverse order (large X first), the integral is carried out from the big end.
  
-  \sa  template <typename Intermediate, typename Final> Final inverse_integrated_density(const std::vector, c2_function &)
- 
-  \param bincenters points at which to sample the c2_function \a binheights
-  \param binheights a c2_function which describes the random number distribution to be produced.
-  \return an interpolating_function  of the type requested in the template which,
-     if evaluated randomly with a uniform variate on [0,1) produces numbers
-     distributed according to \a binheights
-*/
-
-template <typename float_type, typename Final > 
-	Final &inverse_integrated_density(const std::vector<float_type> &bincenters, c2_function<float_type> &binheights)
-{	
-	std::vector<float_type> integral;
-	
-	// integrate from first to last bin in original order, leaving results in integral
-	// ask for relative error of 1e-6 on each bin, with absolute error set to 0 (since we don't know the data scale).
-	float_type sum=binheights.partial_integrals(bincenters, &integral, 0.0, 1e-6); 
-	// the integral vector now has partial integrals... it must be accumulated by summing
-	integral.insert(integral.begin(), 0.0); // integral from start to start is 0
-	float_type scale=1.0/sum;
-	for(size_t i=1; i<integral.size(); i++) integral[i]=integral[i]*scale + integral[i-1];
-	integral.back()=1.0; // force exact value on boundary
-	
-	return  *new Final(integral, bincenters, 
-					   false, 1.0/(scale*binheights(bincenters.front() )), 
-					   false, 1.0/(scale*binheights(bincenters.back() ))
-					   ); // use integral as x axis in inverse function
-}
-
-/**
- \brief Construct a function useful for generation of random numbers from the given distribution
-
- \code
- template <typename Intermediate, typename Final> 
-     Final & inverse_integrated_density(const std::vector &bincenters, const std::vector &binheights) 
- \endcode
- is a variant of \code
- template <typename Final> 
-     Final & inverse_integrated_density(const std::vector &bincenters, c2_function &binheights)
- \endcode
- which takes two std::vectors and generates the intermediate interpolating_function  required for 
- inverse_integrated_density(), and then calls it.
-
- \param bincenters points at which \a binheights are defined
- \param binheights an std::vector which describes the random number distribution to be produced.
- \return an interpolating_function  of the type requested in the template which,
- if evaluated randomly with a uniform variate on [0,1) produces numbers
-  distributed according to \a binheights
-*/
-
-template <typename float_type, typename Intermediate, typename Final> Final
-    &inverse_integrated_density(const std::vector<float_type> &bincenters, const std::vector<float_type> &binheights) 
-{	
-	std::vector<float_type> be(bincenters), bh(binheights);
-	
-	if(be[1] < be[0]) { // reverse data for interpolator if x axis passed in backwards
-		std::reverse(be.begin(), be.end());
-		std::reverse(bh.begin(), bh.end());
-	}
-	
-	Intermediate temp(be, bh); // create a temporary interpolating_function  to integrate
-	Final &result=inverse_integrated_density<Final>(bincenters, temp);
-	
-	return result;
-}
+	\param bincenters the centers of the bins. 
+	\param binheights a c2_function which describes the density of the random number distribution to be produced.
+	\return an interpolating_function_p  of the type requested in the template which,
+	if evaluated randomly with a uniform variate on [0,1] produces numbers
+	distributed according to \a binheights
+ */
+template <typename float_type, typename Final> 
+interpolating_function_p<float_type>  & inverse_integrated_density_function(
+		const std::vector<float_type> &bincenters, const c2_function<float_type> &binheights)
+		throw(c2_exception);
 
 /// \brief create a c2_function which smoothly connects two other c2_functions.
-///
+/// \ingroup parametric_functions
 /// This takes two points and generates a polynomial which matches two c2_function arguments 
 /// at those two points, with two derivatives at each point, and an arbitrary value at the center of the 
@@ -1414,40 +2018,88 @@
 /// of order 5.  If \a auto_center is false, the value \a y1 is used at the midpoint, resulting in a 
 /// polynomial of order 6.
-template <typename float_type=double> class c2_connector_function : public c2_function<float_type> {
+///
+/// This is usually used in conjunction with c2_piecewise_function_p to assemble an apparently seamless 
+/// function from a series of segments. 
+/// \see \ref piecewise_applications_subsec "Sample Applications" and \ref c2_function::adaptively_sample() "Adaptive sampling"
+///
+/// The factory function c2_factory::connector_function() creates *new c2_connector_function_p
+template <typename float_type=double> class c2_connector_function_p : public c2_function<float_type> {
 public:	
-	/// \brief construct the container
-	/// \param f1 the function on the left side to be connected 
+	/// \brief construct the container from two functions
+	/// \param x0 the point at which to match \a f1 and its derivatives
+	/// \param f0 the function on the left side to be connected 
+	/// \param x2 the point at which to match \a f2 and its derivatives
 	/// \param f2 the function on the right side to be connected
-	/// \param x0 the point at which to match \a f1 and its derivatives
-	/// \param x2 the point at which to match \a f2 and its derivatives
 	/// \param auto_center if true, no midpoint value is specified.  If false, match the value \a y1 at the midpoint
 	/// \param y1 the value to match at the midpoint, if \a auto_center is false
-	/// \return a c2_function with domain (\a x0,\a x2) which smoothly connects \a f1 and \a f2
-	c2_connector_function(const c2_function<float_type> &f1, const c2_function<float_type> &f2, float_type x0, float_type x2, 
+	/// \return a c2_function with domain (\a x0,\a x2) which smoothly connects \a f0(x0) and \a f2(x2)
+	c2_connector_function_p(float_type x0, const c2_function<float_type> &f0, float_type x2, const c2_function<float_type> &f2, 
 			bool auto_center, float_type y1);
+	/// \brief construct the container from numerical values
+	/// \param x0 the position of the left edge
+	/// \param y0 the function derivative on the left boundary
+	/// \param yp0 the function second derivative on the left boundary
+	/// \param ypp0 the function value on the left boundary
+	/// \param x2 the position of the right edge
+	/// \param y2 the function derivative on the right boundary
+	/// \param yp2 the function second derivative on the right boundary
+	/// \param ypp2 the function value on the right boundary
+	/// \param auto_center if true, no midpoint value is specified.  If false, match the value \a y1 at the midpoint
+	/// \param y1 the value to match at the midpoint, if \a auto_center is false
+	/// \return a c2_function with domain (\a x0,\a x2) which smoothly connects the points described
+	/// \anchor c2_connector_raw_init_docs
+	c2_connector_function_p(
+		float_type x0, float_type y0, float_type yp0, float_type ypp0, 
+		float_type x2, float_type y2, float_type yp2, float_type ypp2, 
+		bool auto_center, float_type y1);
+	/// \brief construct the container from c2_fblock<float_type> objects
+	/// \param fb0 the left edge
+	/// \param fb2 the right edge
+	/// \param auto_center if true, no midpoint value is specified.  If false, match the value \a y1 at the midpoint
+	/// \param y1 the value to match at the midpoint, if \a auto_center is false
+	/// \return a c2_function with domain (\a fb0.x,\a fb2.x) which smoothly connects \a fb0 and \a fb2
+	c2_connector_function_p(
+		const c2_fblock<float_type> &fb0, 
+		const c2_fblock<float_type> &fb2, 
+		bool auto_center, float_type y1);
+
 	/// \brief destructor
-	virtual ~c2_connector_function();
+	virtual ~c2_connector_function_p();
 	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception);
 protected:
-	float_type fx1, fhinv, fdx, fy1, fa, fb, fc, fd, fe, ff;
-};
-
-
+	/// \brief fill container numerically
+	void init(
+		const c2_fblock<float_type> &fb0, 
+		const c2_fblock<float_type> &fb2, 
+		bool auto_center, float_type y1);
+
+	float_type fhinv, fy1, fa, fb, fc, fd, fe, ff;
+};
 
 /// \brief create a c2_function which is a piecewise assembly of other c2_functions.
-///
+/// \ingroup containers
 /// The functions must have increasing, non-overlapping domains.  Any empty space
 /// between functions will be filled with a linear interpolation.
+///
+/// \note If you want a smooth connection, instead of the default linear interpolation,
+/// create a c2_connector_function_p to bridge the gap.  The linear interpolation is intended 
+/// to be a barely intelligent bridge, and may never get used by anyone.
+///
 /// \note The creation of the container results in the creation of an explicit sampling grid.  
 /// If this is used with functions with a large domain, or which generate very dense sampling grids,
 /// it could eat a lot of memory.  Do not abuse this by using functions which can generate gigantic grids.
 /// 
-/// See c2_plugin_function for a discussion of how this might be used.
-template <typename float_type=double> class c2_piecewise_function : public c2_function<float_type> {
+/// \see \ref piecewise_applications_subsec "Sample Applications" \n
+/// c2_plugin_function_p page \n
+/// c2_connector_function_p page \n
+/// \ref c2_function::adaptively_sample() "Adaptive sampling"
+///
+/// The factory function c2_factory::piecewise_function() creates *new c2_piecewise_function_p
+template <typename float_type=double> class c2_piecewise_function_p : public c2_function<float_type> {
 public:	
 	/// \brief construct the container
-	c2_piecewise_function();
+	c2_piecewise_function_p();
 	/// \brief destructor
-	virtual ~c2_piecewise_function();
+	virtual ~c2_piecewise_function_p();
 	virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception);
 	/// \brief append a new function to the sequence
@@ -1457,11 +2109,8 @@
 	/// the final function.  If the domain exactly abuts the domain of the previous function, it
 	/// will be directly attached.  If there is a gap, the gap will be filled in by linear interpolation.
-	/// If the function being appended is to be deleted automatically when this container is deleted, set the pass_ownership flag.
 	/// \param func a c2_function with a defined domain to be appended to the collection
-	/// \param pass_ownership if set, \a func will be deleted when the container is destroyed
-	void append_function(c2_function<float_type> &func, bool pass_ownership) throw (c2_exception);
+	void append_function(const c2_function<float_type> &func) throw (c2_exception);
 protected:
-	std::vector<c2_function<float_type> *> functions;
-	std::vector<bool> owns;
+	std::vector<c2_const_ptr<float_type> > functions;
 	mutable int lastKLow;
 };
Index: trunk/examples/extended/electromagnetic/TestEm7/include/c2_function.icc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm7/include/c2_function.icc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm7/include/c2_function.icc	(revision 1230)
@@ -1,2 +1,28 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+//
 /**
  *  \file 
@@ -7,5 +33,5 @@
  *  \author Copyright 2005 __Vanderbilt University__. All rights reserved.
  *
- *	\version c2_function.cc,v 1.43 2007/11/12 20:22:54 marcus Exp
+ *	\version c2_function.cc,v 1.169 2008/05/22 12:45:19 marcus Exp
  */
 
@@ -22,5 +48,5 @@
 
 template <typename float_type> const std::string c2_function<float_type>::cvs_file_vers() const 
-{ return "c2_function.cc,v 1.43 2007/11/12 20:22:54 marcus Exp"; }
+{ return "c2_function.cc,v 1.169 2008/05/22 12:45:19 marcus Exp"; }
 
 // find a pre-bracketed root of a c2_function, which is a MUCH easier job than general root finding
@@ -33,5 +59,4 @@
 	// find f(x)=value within the brackets, using the guarantees of smoothness associated with a c2_function
 	// can use local f(x)=a*x**2 + b*x + c and solve quadratic to find root, then iterate
-	reset_evaluations();
 	
 	float_type yp, yp2; // we will make unused pointers point here, to save null checks later
@@ -47,22 +72,36 @@
 	float_type c, b;
 	
+	if(!root_info) {
+		root_info=new struct c2_root_info;
+		root_info->inited=false;
+	}
 	// this new logic is to keep track of where we were before, and lower the number of 
 	// function evaluations if we are searching inside the same bracket as before.
 	// Since this root finder has, very often, the bracket of the entire domain of the function,
 	// this makes a big difference, especially to c2_inverse_function
-	if(!rootInitialized || upper_bracket != lastRootUpperX || lower_bracket != lastRootLowerX) {
-		lastRootUpperY=value_with_derivatives(upper_bracket, final_yprime, final_yprime2);
-		increment_evaluations();
-		lastRootUpperX=upper_bracket;
-	
-		lastRootLowerY=value_with_derivatives(lower_bracket, final_yprime, final_yprime2);
-		increment_evaluations();
-		lastRootLowerX=lower_bracket;
-		rootInitialized=true;
-	}
-	
-	float_type clower=lastRootLowerY-value;
-	float_type cupper=lastRootUpperY-value;
-	if(clower*cupper >0) { 
+	if(!root_info->inited || upper_bracket != root_info->upper.x || lower_bracket != root_info->lower.x) {
+		root_info->upper.x=upper_bracket;
+		fill_fblock(root_info->upper);		
+		root_info->lower.x=lower_bracket;
+		fill_fblock(root_info->lower);
+		root_info->inited=true;
+	}
+	
+	float_type clower=root_info->lower.y-value;
+	if(!clower) {
+		*final_yprime=root_info->lower.yp;
+		*final_yprime2=root_info->lower.ypp;
+		return lower_bracket;
+	}
+	
+	float_type cupper=root_info->upper.y-value;
+	if(!cupper) {
+		*final_yprime=root_info->upper.yp;
+		*final_yprime2=root_info->upper.ypp;
+		return upper_bracket;
+	}
+	const float_type lower_sign = (clower < 0) ? -1 : 1;
+	
+	if(lower_sign*cupper >0) { 
 		// argh, no sign change in here!
 		if(error) { *error=1; return 0.0; }
@@ -104,9 +143,13 @@
 	   }
 	   c=value_with_derivatives(root, final_yprime, final_yprime2)-value; // compute initial values
+	   if(c2_isnan(c)) {
+		   bad_x_point=root;
+		   return c; // return the nan if a computation failed
+	   }
 	   b=*final_yprime; // make a local copy for readability
 	   increment_evaluations();
 
 	   // now, close in bracket on whichever side this still brackets
-	   if(c*clower < 0.0) {
+	   if(c*lower_sign < 0.0) {
 		   cupper=c;
 		   upper_bracket=root;
@@ -131,18 +174,22 @@
 
 // the recursive part of the integrator is agressively designed to minimize copying of data... lots of pointers
-template <typename float_type> float_type c2_function<float_type>::integrate_step(c2_integrate_recur &rb) const
-{
-	struct c2_integrate_fblock *fbl[3]={rb.f0, rb.f1, rb.f2}; 
-	struct c2_integrate_fblock f1; // will hold new middle values
-	float_type retvals[2]={0.0,0.0};
-	float_type lr[2];
-
-	// std::cout << "entering with " << rb.f0->x << " " << rb.f1->x << " " << rb.f2->x << std::endl;
-	
-	int depth=rb.depth; // save this from the recursion block
-	float_type abs_tol=rb.abs_tol;  // this is the value we will pass down
-	float_type *rblr=rb.lr; // save pointer to our parent's lr[2] array since it will get trampled in recursion
-	
-	if(!depth) {
+template <typename float_type> float_type c2_function<float_type>::integrate_step(c2_integrate_recur &rb) const throw(c2_exception)
+{
+	std::vector< recur_item > &rb_stack=*rb.rb_stack; // heap-based stack of data for recursion
+	rb_stack.clear();
+
+	recur_item top;
+	top.depth=0; top.done=false; top.f0index=0; top.f2index=0; top.step_sum=0;
+	
+	// push storage for our initial elements
+	rb_stack.push_back(top);
+	rb_stack.back().f1=*rb.f0;
+	rb_stack.back().done=true; // this element will never be evaluated further
+	
+	rb_stack.push_back(top);
+	rb_stack.back().f1=*rb.f1;
+	rb_stack.back().done=true; // this element will never be evaluated further
+		
+	if(!rb.inited) {
 		switch(rb.derivs) {
 			case 0:
@@ -157,57 +204,97 @@
 		
 		rb.extrap2=1.0/(rb.extrap_coef-1.0);
-	}
-	
-	for (int i=0; i<(depth==0?1:2); i++) { // handle left and right intervals, but only left one for depth=0
-		struct c2_integrate_fblock *f0=fbl[i], *f2=fbl[i+1];
-		f1.x=0.5*(f0->x + f2->x); // center of interval
-		float_type dx=f2->x - f0->x;
-		float_type dx2 = 0.5*dx;
-		float_type total;
-		
-		f1.y=value_with_derivatives(f1.x, &(f1.yp), &(f1.ypp));
-		increment_evaluations(); 
-		
+		rb.dx_tolerance=10.0*std::numeric_limits<float_type>::epsilon();
+		rb.abs_tol_min=10.0*std::numeric_limits<float_type>::min();
+		rb.inited=true;
+	}
+	
+	// now, push our first real element
+	top.f0index=0; // left element is stack[0]
+	top.f2index=1; // right element is stack[1]
+	top.abs_tol=rb.abs_tol;
+	rb_stack.push_back(top);
+		
+	while(rb_stack.size() > 2) {
+		recur_item &back=rb_stack.back();
+		if(back.done) {
+			float_type sum=back.step_sum;
+			rb_stack.pop_back();
+			rb_stack.back().step_sum+=sum; // bump our sum up to the parent
+			continue;
+		}
+		back.done=true;
+		
+		c2_fblock<float_type> &f0=rb_stack[back.f0index].f1, &f2=rb_stack[back.f2index].f1; 
+		c2_fblock<float_type> &f1=back.f1; // will hold new middle values
+		size_t f1index=rb_stack.size()-1; // our current offset
+		float_type abs_tol=back.abs_tol;
+		
+		f1.x=0.5*(f0.x + f2.x); // center of interval
+		float_type dx2=0.5*(f2.x - f0.x);
+
 		// check for underflow on step size, which prevents us from achieving specified accuracy. 
-		if(std::abs(dx) < std::abs(f1.x)*rb.rel_tol) {
+		if(std::abs(dx2) < std::abs(f1.x)*rb.dx_tolerance || std::abs(dx2) < rb.abs_tol_min) {
 			std::ostringstream outstr;
-			outstr << "Step size underflow in adaptive_partial_integrals at depth=" << depth << ", x= " << f1.x;
+			outstr << "Step size underflow in adaptive_partial_integrals at depth=" << back.depth << ", x= " << f1.x;
 			throw c2_exception(outstr.str().c_str());
 		}
 		
-		if(!depth) { // top level, total has not been initialized yet
-			switch(rb.derivs) { // create estimate of next lower order for first try
+		fill_fblock(f1);
+		if(c2_isnan(f1.y)) {
+			bad_x_point=f1.x;
+			return f1.y; // can't go any further if a nan has appeared
+		}
+		
+		bool yptrouble=f0.ypbad || f2.ypbad || f1.ypbad;
+		bool ypptrouble=f0.yppbad || f2.yppbad || f1.yppbad;
+		
+		// select the real derivative count based on whether we are at a point where derivatives exist
+		int derivs = std::min(rb.derivs, (yptrouble||ypptrouble)?(yptrouble?0:1):2);
+				
+		if(!back.depth) { // top level, total has not been initialized yet
+			switch(derivs) { // create estimate of next lower order for first try
 				case 0:
-					total=0.5*(f0->y+f2->y)*dx; break;
+					back.previous_estimate=(f0.y+f2.y)*dx2; break;
 				case 1:
-					total=(f0->y+4.0*f1.y+f2->y)*dx/6.0; break;
+					back.previous_estimate=(f0.y+4.0*f1.y+f2.y)*dx2/3.0; break;
 				case 2:
-					total=( (14*f0->y + 32*f1.y + 14*f2->y) +  dx * (f0->yp - f2->yp) ) * dx /60.; break;
+					back.previous_estimate=( (14*f0.y + 32*f1.y + 14*f2.y) +  2*dx2 * (f0.yp - f2.yp) ) * dx2 /30.; break;
 				default:
-					total=0.0; // just to suppress missing default warnings
+					back.previous_estimate=0.0; // just to suppress missing default warnings
 			}
-		} else total=rblr[i]; // otherwise, get it from previous level
+		} 
 		
 		float_type left, right;
 		
-		switch(rb.derivs) {
+		// pre-compute constants so all multiplies use a small dynamic range
+		// constants for 0 derivative integrator
+		static const float_type c0c1=5./12., c0c2=8./12., c0c3=-1./12.;
+		// constants for 1 derivative integrator
+		static const float_type c1c1=101./240., c1c2=128./240., c1c3=11./240.,
+			c1c4=13./240., c1c5=-40./240., c1c6=-3./240.;
+		// constants for 2 derivative integrator
+		static const float_type c2c1=169./40320., c2c2=1024./ 40320., c2c3=-41./40320.,
+			c2c4=2727./40320., c2c5=-5040./40320., c2c6=423./40320.,
+			c2c7=17007./40320., c2c8=24576./40320., c2c9=-1263./40320.;
+		
+		switch(derivs) {
 			case 2:
 				// use ninth-order estimates for each side, from full set of all values (!) (Thanks, Mathematica!)
-				left=	( ( (169*f0->ypp + 1024*f1.ypp -  41*f2->ypp)*dx2 + 
-							(2727*f0->yp - 5040*f1.yp +  423*f2->yp) )*dx2 + 
-						  (17007*f0->y + 24576*f1.y -  1263*f2->y) )* (dx2/40320.0);
-				right=	( ( (169*f2->ypp + 1024*f1.ypp -  41*f0->ypp)*dx2 - 
-							(2727*f2->yp - 5040*f1.yp +  423*f0->yp) )*dx2 + 
-						  (17007*f2->y + 24576*f1.y -  1263*f0->y) )* (dx2/40320.0);
-				// std::cout << f0->x << " " << f1.x << " " << f2->x <<  std::endl ;
-				// std::cout << f0->y << " " << f1.y << " " << f2->y << " " << left << " " << right << " " << total << std::endl ;
+				left=	( ( (c2c1*f0.ypp + c2c2*f1.ypp +  c2c3*f2.ypp)*dx2 + 
+							(c2c4*f0.yp + c2c5*f1.yp +  c2c6*f2.yp) )*dx2 + 
+						  (c2c7*f0.y + c2c8*f1.y +  c2c9*f2.y) )* dx2;
+				right=	( ( (c2c1*f2.ypp + c2c2*f1.ypp +  c2c3*f0.ypp)*dx2 - 
+							(c2c4*f2.yp + c2c5*f1.yp +  c2c6*f0.yp) )*dx2 + 
+						  (c2c7*f2.y + c2c8*f1.y +  c2c9*f0.y) )* dx2;
+				// std::cout << f0.x << " " << f1.x << " " << f2.x <<  std::endl ;
+				// std::cout << f0.y << " " << f1.y << " " << f2.y << " " << left << " " << right << " " << total << std::endl ;
 				break;
 			case 1:
-				left= 	( (202*f0->y + 256*f1.y + 22*f2->y) + dx*(13*f0->yp - 40*f1.yp - 3*f2->yp) ) * dx /960.;
-				right= 	( (202*f2->y + 256*f1.y + 22*f0->y) - dx*(13*f2->yp - 40*f1.yp - 3*f0->yp) ) * dx /960.;
+				left= 	( (c1c1*f0.y + c1c2*f1.y + c1c3*f2.y) + dx2*(c1c4*f0.yp + c1c5*f1.yp + c1c6*f2.yp) ) * dx2 ;
+				right= 	( (c1c1*f2.y + c1c2*f1.y + c1c3*f0.y) - dx2*(c1c4*f2.yp + c1c5*f1.yp + c1c6*f0.yp) ) * dx2 ;
 				break;
 			case 0:
-				left=	(5*f0->y + 8*f1.y - f2->y)*dx/24.;
-				right=	(5*f2->y + 8*f1.y - f0->y)*dx/24.;
+				left=	(c0c1*f0.y + c0c2*f1.y + c0c3*f2.y)*dx2;
+				right=	(c0c1*f2.y + c0c2*f1.y + c0c3*f0.y)*dx2;
 				break;
 			default:
@@ -216,32 +303,40 @@
 		}
 		
-		lr[0]= left; // left interval
-		lr[1]= right; // right interval
 		float_type lrsum=left+right;
 
-		float_type eps=std::abs(total-lrsum)*rb.eps_scale; 
-		if(rb.extrapolate) eps*=rb.eps_scale; 
-		
-		if(!rb.adapt || eps < abs_tol ||  eps < std::abs(total)*rb.rel_tol) {
-			if(depth==0 || !rb.extrapolate) retvals[i]=lrsum;
-			else {
-				retvals[i]=(rb.extrap_coef*lrsum - total)*rb.extrap2; 	
-				// std::cout << "extrapolating " << lrsum << " " << total << " " << retvals[i] << std::endl;
-				
-			}
+		bool extrapolate=back.depth && rb.extrapolate && (derivs==rb.derivs); // only extrapolate if no trouble with derivs
+		float_type eps=std::abs(back.previous_estimate-lrsum)*rb.eps_scale; 
+		if(extrapolate) eps*=rb.eps_scale; 
+		
+		if(rb.adapt && eps > abs_tol &&  eps > std::abs(lrsum)*rb.rel_tol) {
+			// tolerance not met, subdivide & recur
+			if(abs_tol > rb.abs_tol_min) abs_tol=abs_tol*0.5; // each half has half the error budget
+			top.abs_tol=abs_tol;
+			top.depth=back.depth+1;
+			
+			// save the last things we need from back before a push happens, in case 
+			// the push causes a reallocation and moves the whole stack.
+			size_t f0index=back.f0index, f2index=back.f2index;
+			
+			top.f0index=f1index; top.f2index=f2index; // insert pointers to right side data into our recursion block
+			top.previous_estimate=right;
+			rb_stack.push_back(top);
+
+			top.f0index=f0index; top.f2index=f1index; // insert pointers to left side data into our recursion block
+			top.previous_estimate=left;
+			rb_stack.push_back(top);
+
+		} else if(extrapolate) {
+			// extrapolation only happens on leaf nodes, where the tolerance was met.
+			back.step_sum+=(rb.extrap_coef*lrsum - back.previous_estimate)*rb.extrap2; 
 		} else {
-			rb.depth=depth+1; // increment depth counter
-			rb.lr=lr;  // point to our left-right values array for recursion
-			rb.abs_tol=abs_tol*0.5; // each half has half the error budget
-			rb.f0=f0; rb.f1=&f1; rb.f2=f2; // insert pointers to data into our recursion block
-			// std::cout << "recurring with " << f0->x << " " << f1.x << " " << f2->x <<  std::endl ;
-			retvals[i]=integrate_step(rb); // and recur
-		}
-	}	
-	return retvals[0]+retvals[1];
+			back.step_sum+=lrsum; 
+		}
+	}
+	return rb_stack.back().step_sum; // last element on the stack holds the sum
 }
 
 template <typename float_type> bool c2_function<float_type>::check_monotonicity(
-	const std::vector<float_type> &data, const char message[]) throw(c2_exception)
+	const std::vector<float_type> &data, const char message[]) const throw(c2_exception)
 {
 	size_t np=data.size();
@@ -269,7 +364,9 @@
 }
 
-template <typename float_type> std::vector<float_type> &c2_function<float_type>::get_sampling_grid(float_type xmin, float_type xmax) const 
-{
-	std::vector<float_type> *result=new std::vector<float_type>;
+template <typename float_type> void c2_function<float_type>::
+	get_sampling_grid(float_type xmin, float_type xmax, std::vector<float_type> &grid) const 
+{
+	std::vector<float_type> *result=&grid;
+	result->clear();
 	
 	if( !(sampling_grid) || !(sampling_grid->size()) || (xmax <= sampling_grid->front()) || (xmin >= sampling_grid->back()) ) { 
@@ -313,9 +410,8 @@
 		std::copy(sg.begin()+firstindex, sg.begin()+lastindex+1, result->begin()+initsize);
 		result->back()=xmax;
-		
+
 		//  this is the unrefined sampling grid... now check for very close points on front & back and fix if needed.
 		preen_sampling_grid(result);
 	}
-	return *result;
 }
 
@@ -351,6 +447,6 @@
 }
 
-template <typename float_type> std::vector<float_type> &c2_function<float_type>::
-	refine_sampling_grid(const std::vector<float_type> &grid, size_t refinement) const
+template <typename float_type> void c2_function<float_type>::
+	refine_sampling_grid(std::vector<float_type> &grid, size_t refinement) const
 {
 	size_t np=grid.size();
@@ -358,61 +454,68 @@
 	float_type dxscale=1.0/refinement;
 	
-	std::vector<float_type> *result=new std::vector<float_type>(count);
+	std::vector<float_type> result(count);
 	
 	for(size_t i=0; i<(np-1); i++) {
 		float_type x=grid[i];
 		float_type dx=(grid[i+1]-x)*dxscale;
-		for(size_t j=0; j<refinement; j++, x+=dx) (*result)[i*refinement+j]=x;
-	}
-	(*result)[count-1]=grid.back();
-	return *result;
+		for(size_t j=0; j<refinement; j++, x+=dx) result[i*refinement+j]=x;
+	}
+	result.back()=grid.back();
+	grid=result; // copy the expanded grid back to the input
 }
 
 template <typename float_type> float_type c2_function<float_type>::integral(float_type xmin, float_type xmax, std::vector<float_type> *partials,
-	 float_type abs_tol, float_type rel_tol, int derivs, bool adapt, bool extrapolate) const
-{
-	std::vector<float_type> &grid=get_sampling_grid(xmin, xmax);
-	float_type intg=partial_integrals(grid, partials, abs_tol, rel_tol, adapt, extrapolate);
-	delete &grid;
+	 float_type abs_tol, float_type rel_tol, int derivs, bool adapt, bool extrapolate) const throw(c2_exception)
+{
+	if(xmin==xmax) {
+		if(partials) partials->clear();
+		return 0.0;
+	}
+	std::vector<float_type> grid;
+	get_sampling_grid(xmin, xmax, grid);
+	float_type intg=partial_integrals(grid, partials, abs_tol, rel_tol, derivs, adapt, extrapolate);
 	return intg;
 }
 
 template <typename float_type> c2_function<float_type> &c2_function<float_type>::normalized_function(float_type xmin, float_type xmax, float_type norm) 
+	const throw(c2_exception)
 {
 	float_type intg=integral(xmin, xmax);
-	return *new c2_scaled_function<float_type>(*this, norm/intg);
-}
-
-template <typename float_type> c2_function<float_type> &c2_function<float_type>::square_normalized_function(float_type xmin, float_type xmax, float_type norm) 
-{
-	c2_quadratic<float_type> q(0., 0., 0., 1.);
-	c2_composed_function<float_type> mesquared(q,*this);
-	
-	std::vector<float_type> grid(get_sampling_grid(xmin, xmax));		
-	float_type intg=mesquared.partial_integrals(grid);
-	
-	return *new c2_scaled_function<float_type>(*this, std::sqrt(norm/intg));
+	return *new c2_scaled_function_p<float_type>(*this, norm/intg);
+}
+
+template <typename float_type> c2_function<float_type> &c2_function<float_type>::square_normalized_function(float_type xmin, float_type xmax, float_type norm)
+	const throw(c2_exception)
+{
+	c2_ptr<float_type> mesquared((*new c2_quadratic_p<float_type>(0., 0., 0., 1.))(*this));
+	
+	std::vector<float_type> grid;
+	get_sampling_grid(xmin, xmax, grid);		
+	float_type intg=mesquared->partial_integrals(grid);
+	
+	return *new c2_scaled_function_p<float_type>(*this, std::sqrt(norm/intg));
 }
 
 template <typename float_type> c2_function<float_type> &c2_function<float_type>::square_normalized_function(
 		float_type xmin, float_type xmax, const c2_function<float_type> &weight, float_type norm) 
-{
-	c2_quadratic<float_type> q(0., 0., 0., 1.);
-	c2_composed_function<float_type> mesquared(q,*this);
-	c2_product<float_type> weighted(mesquared, weight); 
-
-	std::vector<float_type> grid(get_sampling_grid(xmin, xmax));	
-	float_type intg=weighted.partial_integrals(grid);
-
-	return *new c2_scaled_function<float_type>(*this, std::sqrt(norm/intg));
+	const throw(c2_exception)
+{
+	c2_ptr<float_type> weighted((*new c2_quadratic_p<float_type>(0., 0., 0., 1.))(*this) * weight);
+
+	std::vector<float_type> grid;
+	get_sampling_grid(xmin, xmax, grid);	
+	float_type intg=weighted->partial_integrals(grid);
+
+	return *new c2_scaled_function_p<float_type>(*this, std::sqrt(norm/intg));
 }
 
 template <typename float_type> float_type c2_function<float_type>::partial_integrals(
 	std::vector<float_type> xgrid, std::vector<float_type> *partials,
-	float_type abs_tol, float_type rel_tol, int derivs, bool adapt, bool extrapolate) const
+	float_type abs_tol, float_type rel_tol, int derivs, bool adapt, bool extrapolate) 
+	const throw(c2_exception)
 {
 	int np=xgrid.size();
 	
-	struct c2_integrate_fblock f0, f2;
+	c2_fblock<float_type> f0, f2;
 	struct c2_integrate_recur rb;
 	rb.rel_tol=rel_tol;
@@ -420,6 +523,9 @@
 	rb.adapt=adapt;
 	rb.derivs=derivs;
-	
-	reset_evaluations(); // counter returns with total evaluations needed for this integral
+	std::vector< recur_item > rb_stack;
+	rb_stack.reserve(20); // enough for most operations
+	rb.rb_stack=&rb_stack;
+	rb.inited=false;
+	float_type dx_inv=1.0/std::abs(xgrid.back()-xgrid.front());
 	
 	if(partials) partials->resize(np-1);
@@ -428,6 +534,9 @@
 	
 	f2.x=xgrid[0];
-	f2.y=value_with_derivatives(f2.x, &f2.yp, &f2.ypp);
-	increment_evaluations();
+	fill_fblock(f2);
+	if(c2_isnan(f2.y)) {
+		bad_x_point=f2.x;
+		return f2.y; // can't go any further if a nan has appeared
+	}
 	
 	for(int i=0; i<np-1; i++) {
@@ -435,63 +544,69 @@
 		
 		f2.x=xgrid[i+1];
-		f2.y=value_with_derivatives(f2.x, &f2.yp, &f2.ypp);
-		increment_evaluations();
-		
-		rb.depth=0;
-		rb.abs_tol=abs_tol;
-		rb.f0=&f0; rb.f1=&f2; rb.f2=&f2; // we are really only using the left half for the top level
-		rb.lr=0; // pointer is meaningless; will be filled in in recursion
+		fill_fblock(f2);
+		if(c2_isnan(f2.y)) {
+			bad_x_point=f2.x;
+			return f2.y; // can't go any further if a nan has appeared
+		}
+		
+		rb.abs_tol=abs_tol*std::abs(f2.x-f0.x)*dx_inv; // distribute error tolerance over whole domain
+		rb.f0=&f0; rb.f1=&f2; 
 		float_type ps=integrate_step(rb);
 		sum+=ps;
 		if(partials) (*partials)[i]=ps;
+		if(c2_isnan(ps)) break; // NaN stops integration
 	}
 	return sum;
 }
-
-// declare singleton functions for most common c2_function instances
-#define c2_singleton(X) template <typename float_type> const c2_##X<float_type> c2_##X<float_type>::X=c2_##X();
-c2_singleton(sin)
-c2_singleton(cos)
-c2_singleton(tan)
-c2_singleton(log)
-c2_singleton(exp)
-c2_singleton(sqrt)
-c2_singleton(identity)
-
-// reciprocal is actually parametric (a/x), but make singleton 1/x
-template <typename float_type> const c2_recip<float_type> c2_recip<float_type>::recip=c2_recip(1.0);
-
-#undef c2_singleton
 
 // generate a sampling grid at points separated by dx=5, which is intentionally
 // incommensurate with pi and 2*pi so grid errors are somewhat randomized
-template <typename float_type> std::vector<float_type> &c2_sin<float_type>::get_sampling_grid(float_type xmin, float_type xmax)
-{
-	std::vector<float_type> *result=new std::vector<float_type>;
-	
-	for(; xmin < xmax; xmin+=5.0) result->push_back(xmin);
-	result->push_back(xmax);
-	this->preen_sampling_grid(result);
-	return *result;
-}
-
-template <typename float_type> float_type Identity(float_type x) { return x; } // a useful function
-template <typename float_type> float_type f_one(float_type) { return 1.0; } // the first derivative of identity
-template <typename float_type> float_type f_zero(float_type) { return 0.0; } // the second derivative of identity
+template <typename float_type> void c2_sin_p<float_type>::
+	get_sampling_grid(float_type xmin, float_type xmax,  std::vector<float_type> &grid) const
+{
+	grid.clear();	
+	for(; xmin < xmax; xmin+=5.0) grid.push_back(xmin);
+	grid.push_back(xmax);
+	this->preen_sampling_grid(&grid);
+}
+
+template <typename float_type> float_type c2_function_transformation<float_type>::evaluate(
+		float_type xraw, 
+		float_type y, float_type yp0, float_type ypp0,
+		float_type *yprime, float_type *yprime2) const
+{
+	y=Y.fHasStaticTransforms ? Y.pOut(y) : Y.fOut(y); 
+		
+	if(yprime || yprime2) {
+
+		float_type yp, yp2;
+		if(X.fHasStaticTransforms && Y.fHasStaticTransforms) {
+			float_type fpi=1.0/Y.pInPrime(y);
+			float_type gp=X.pInPrime(xraw);
+			// from Mathematica Dt[InverseFunction[f][y[g[x]]], x]
+			yp=gp*yp0*fpi; // transformed derivative
+			yp2=(gp*gp*ypp0 + X.pInDPrime(xraw)*yp0  - Y.pInDPrime(y)*yp*yp )*fpi; 
+		} else {
+			float_type fpi=1.0/Y.fInPrime(y);
+			float_type gp=X.fInPrime(xraw);
+			// from Mathematica Dt[InverseFunction[f][y[g[x]]], x]
+			yp=gp*yp0*fpi; // transformed derivative
+			yp2=(gp*gp*ypp0 + X.fInDPrime(xraw)*yp0  - Y.fInDPrime(y)*yp*yp )*fpi; 
+		} 
+		if(yprime) *yprime=yp;
+		if(yprime2) *yprime2=yp2;
+	}	
+	return y;
+}
 
 //  The constructor
-template <typename float_type> void interpolating_function<float_type>::init(
+template <typename float_type> interpolating_function_p<float_type> & interpolating_function_p<float_type>::load(
 								const std::vector<float_type> &x, const std::vector<float_type> &f, 
 								bool lowerSlopeNatural, float_type lowerSlope, 
 								bool upperSlopeNatural, float_type upperSlope,
-								float_type (*inputXConversion)(float_type), 
-								float_type (*inputXConversionPrime)(float_type), 
-								float_type (*inputXConversionDPrime)(float_type), 
-								float_type (*inputYConversion)(float_type), 
-								float_type (*inputYConversionPrime)(float_type), 
-								float_type (*inputYConversionDPrime)(float_type), 
-								float_type (*outputYConversion)(float_type)
-					) throw(c2_exception) 
-{
+								bool splined
+		) throw(c2_exception) 
+{
+	c2_ptr<float_type> keepme(*this);
 	X= x;
 	F= f;
@@ -501,13 +616,5 @@
 
 	set_domain(std::min(Xraw.front(), Xraw.back()),std::max(Xraw.front(), Xraw.back()));
-	
-	fXin=inputXConversion;
-	fXinPrime=inputXConversionPrime;
-	fXinDPrime=inputXConversionDPrime;
-	fYin=inputYConversion;
-	fYinPrime=inputYConversionPrime;
-	fYinDPrime=inputYConversionDPrime;
-	fYout=outputYConversion;
-	
+		
 	if(x.size() != f.size()) {
 		throw c2_exception("interpolating_function::init() -- x & y inputs are of different size");
@@ -529,23 +636,13 @@
 	}
 	
-	if(fXin) { // check if X scale is nonlinear, and if so, do transform
-		if(!lowerSlopeNatural) lowerSlope /= fXinPrime(X[0]);
-		if(!upperSlopeNatural) upperSlope /= fXinPrime(X[np-1]);
-		for(size_t i=0; i<np; i++) X[i]=fXin(X[i]);
-	} else {
-		fXin=Identity<float_type>;
-		fXinPrime=f_one<float_type>;
-		fXinDPrime=f_zero<float_type>;
-	}
-	
-	if(inputYConversion) {  // check if Y scale is nonlinear, and if so, do transform
-		if(!lowerSlopeNatural) lowerSlope *= fYinPrime(F[0]);
-		if(!upperSlopeNatural) upperSlope *= fYinPrime(F[np-1]);
-		for(size_t i=0; i<np; i++) F[i]=inputYConversion(F[i]);
-	} else {
-		fYin=Identity<float_type>;
-		fYinPrime=f_one<float_type>;
-		fYinDPrime=f_zero<float_type>;
-		fYout=Identity<float_type>;
+	if(fTransform.X.fTransformed) { // check if X scale is nonlinear, and if so, do transform
+		if(!lowerSlopeNatural) lowerSlope /= fTransform.X.fInPrime(X[0]);
+		if(!upperSlopeNatural) upperSlope /= fTransform.X.fInPrime(X[np-1]);
+		for(size_t i=0; i<np; i++) X[i]=fTransform.X.fIn(X[i]);
+	}	
+	if(fTransform.Y.fTransformed) {  // check if Y scale is nonlinear, and if so, do transform
+		if(!lowerSlopeNatural) lowerSlope *= fTransform.Y.fInPrime(F[0]);
+		if(!upperSlopeNatural) upperSlope *= fTransform.Y.fInPrime(F[np-1]);
+		for(size_t i=0; i<np; i++) F[i]=fTransform.Y.fIn(F[i]);
 	}
 		  
@@ -553,7 +650,126 @@
 		"interpolating_function::init() non-monotonic transformed x input"); 
 	
-	// construct spline tables here.  
+	if(splined) spline(lowerSlopeNatural, lowerSlope, upperSlopeNatural, upperSlope);
+	else y2.assign(np,0.0);
+	
+	lastKLow=0; 
+	keepme.release_for_return();
+	return *this;
+}
+
+/*
+//  The constructor
+template <typename float_type> interpolating_function_p<float_type> & interpolating_function_p<float_type>::load_pairs(
+								std::vector<std::pair<float_type, float_type> > &data, 
+								bool lowerSlopeNatural, float_type lowerSlope, 
+								bool upperSlopeNatural, float_type upperSlope,
+								bool splined
+		) throw(c2_exception) 
+{
+	c2_ptr<float_type> keepme(*this);
+	
+	size_t np=data.size();
+	if(np < 2) {
+		throw c2_exception("interpolating_function::init() -- input < 2 elements ");
+	}
+	
+	// sort into ascending order
+	std::sort(data.begin(), data.end(), comp_pair);	
+	
+	std::vector<float_type> xtmp, ytmp;
+	xtmp.reserve(np);
+	ytmp.reserve(np);
+	for (size_t i=0; i<np; i++) {
+		xtmp.push_back(data[i].first);
+		ytmp.push_back(data[i].second);
+	}
+	this->load(xtmp, ytmp, lowerSlopeNatural, lowerSlope, upperSlopeNatural, upperSlope, splined);
+	
+	keepme.release_for_return();
+	return *this;
+}
+
+template <typename float_type> interpolating_function_p<float_type> & 
+	interpolating_function_p<float_type>::load_random_generator_function(
+	const std::vector<float_type> &bincenters, const c2_function<float_type> &binheights)
+	throw(c2_exception)
+{	
+	c2_ptr<float_type> keepme(*this);
+
+	std::vector<float_type> integral;
+	c2_const_ptr<float_type> keepit(binheights); // manage function... not really needed here, but always safe.
+	// integrate from first to last bin in original order, leaving results in integral
+	// ask for relative error of 1e-6 on each bin, with absolute error set to 0 (since we don't know the data scale).
+	float_type sum=binheights.partial_integrals(bincenters, &integral, 0.0, 1e-6); 
+	// the integral vector now has partial integrals... it must be accumulated by summing
+	integral.insert(integral.begin(), 0.0); // integral from start to start is 0
+	float_type scale=1.0/sum;
+	for(size_t i=1; i<integral.size(); i++) integral[i]=integral[i]*scale + integral[i-1];
+	integral.back()=1.0; // force exact value on boundary
+	
+	this->load(integral, bincenters, 
+					   false, 1.0/(scale*binheights(bincenters.front() )), 
+					   false, 1.0/(scale*binheights(bincenters.back() ))
+					   ); // use integral as x axis in inverse function
+	keepme.release_for_return();
+	return *this;
+}
+
+template <typename float_type> interpolating_function_p<float_type> & 
+	interpolating_function_p<float_type>::load_random_generator_bins(
+	const std::vector<float_type> &bins, const std::vector<float_type> &binheights)
+	throw(c2_exception)
+{	
+	c2_ptr<float_type> keepme(*this);
+
+	size_t np=binheights.size();
+	std::vector<float_type> integral(np+1), bin_edges(np+1);
+	
+	// compute the integral based on estimates of the bin edges from the given bin centers...
+	// except for bin 0 & final bin, the edge of a bin is halfway between then center of the 
+	// bin and the center of the previous/next bin.
+	// This gives width[n] = (center[n+1]+center[n])/2 - (center[n]+center[n-1])/2 = (center[n+1]-center[n-1])/2
+	// for the edges, assume a bin of width (center[1]-center[0]) or (center[np-1]-center[np-2])
+	// be careful that absolute values are used in case data are reversed.
+
+	if(bins.size() == binheights.size()+1) {
+		bin_edges=bins; // edges array was passed in
+	} else if (bins.size() == binheights.size()) {
+		bin_edges.front()=bins[0] - (bins[1]-bins[0])*0.5; // edge bin
+		for(size_t i=1; i<np; i++) {
+			bin_edges[i]=(bins[i]+bins[i-1])*0.5;
+		}
+		bin_edges.back()=bins[np-1] + (bins[np-1]-bins[np-2])*0.5; // edge bin
+	} else {
+		throw c2_exception("inconsistent bin vectors passed to load_random_generator_bins");
+	}
+	
+	float_type running_sum=0.0;
+	for(size_t i=0; i<np; i++) {
+		integral[i]=running_sum;
+		if(!binheights[i]) throw c2_exception("empty bin passed to load_random_generator_bins");
+		running_sum+=binheights[i]*std::abs(bin_edges[i+1]-bin_edges[i]);
+	}
+	float_type scale=1.0/running_sum;
+	for(size_t i=0; i<np; i++) integral[i]*=scale;
+	integral.back()=1.0; // force exactly correct value on boundary
+	this->load(integral, bin_edges, 
+		  false, 1.0/(scale*binheights.front()), 
+		  false, 1.0/(scale*binheights.back())
+		  ); // use integral as x axis in inverse function
+	keepme.release_for_return();
+	return *this;
+}
+*/
+
+//  The spline table generator
+template <typename float_type> void interpolating_function_p<float_type>::spline(
+	 bool lowerSlopeNatural, float_type lowerSlope, 
+	 bool upperSlopeNatural, float_type upperSlope
+	 ) throw(c2_exception) 
+{
+// construct spline tables here.  
 	// this code is a re-translation of the pythonlabtools spline algorithm from pythonlabtools.sourceforge.net
-	
+	size_t np=X.size();
 	std::vector<float_type> u(np),  dy(np-1), dx(np-1), dxi(np-1), dx2i(np-2), siga(np-2), dydx(np-1);
 	
@@ -596,13 +812,102 @@
 	y2[np-1]=(un-qn*u[np-2])/(qn*y2[np-2]+1.0);
 	for (size_t k=np-1; k != 0; k--) y2[k-1]=y2[k-1]*y2[k]+u[k-1];
-	
-	lastKLow=-1; // flag  new X search required for next evaluation
+}
+
+template <typename float_type> interpolating_function_p<float_type> &interpolating_function_p<float_type>::sample_function(
+			const c2_function<float_type> &func,
+			float_type xmin, float_type xmax, float_type abs_tol, float_type rel_tol, 
+			bool lowerSlopeNatural, float_type lowerSlope, 
+			bool upperSlopeNatural, float_type upperSlope
+ ) throw(c2_exception) 
+{ 
+	c2_ptr<float_type> keepme(*this);
+	
+	const c2_transformation<float_type> &XX=fTransform.X, &YY=fTransform.Y; // shortcuts
+	
+	 // set up our params to look like the samplng function for now
+	 sampler_function=func;
+	 std::vector<float_type> grid;
+	 func.get_sampling_grid(xmin, xmax, grid);
+	 size_t gsize=grid.size();
+	 if(XX.fTransformed) for(size_t i=0; i<gsize; i++) grid[i]=XX.fIn(grid[i]);
+	 set_sampling_grid_pointer(grid);
+	 
+	// float_type xmin1=fXin(xmin), xmax1=fXin(xmax); // bounds in transformed space
+	 // get a list of points needed in transformed space, directly into our tables
+	 this->adaptively_sample(grid.front(), grid.back(), 8*abs_tol, 8*rel_tol, 0, &X, &F);
+	 // clear the sampler function now, since otherwise our value_with_derivatives is broken
+	 sampler_function.unset_function();
+	 
+	 xInverted=check_monotonicity(X, 
+		  "interpolating_function::init() non-monotonic transformed x input"); 
+	 
+	 size_t np=X.size();
+	 
+	 // Xraw is useful in some of the arithmetic operations between interpolating functions
+	 if(!XX.fTransformed) Xraw=X;
+	 else {
+		 Xraw.resize(np);
+		 for (size_t i=1; i<np-1; i++) Xraw[i]=XX.fOut(X[i]);
+		 Xraw.front()=xmin;
+		 Xraw.back()=xmax;
+	 }
+	 	 
+	 bool xraw_rev=check_monotonicity(Xraw, 
+		  "interpolating_function::init() non-monotonic raw x input"); 
+		  // which way does raw X point?  sampling grid MUST be increasing
+	 
+	 if(!xraw_rev) { // we can use pointer to raw X values if they are in the right order
+		 set_sampling_grid_pointer(Xraw); 
+		 // our intial grid of x values is certainly a good guess for 'interesting' points
+	 } else {
+		 set_sampling_grid(Xraw); // make a copy of it, and assure it is in right order
+	 }
+	 
+	 if(XX.fTransformed) { // check if X scale is nonlinear, and if so, do transform
+		 if(!lowerSlopeNatural) lowerSlope /= XX.fInPrime(xmin);
+		 if(!upperSlopeNatural) upperSlope /= XX.fInPrime(xmax);
+	 }	
+	 if(YY.fTransformed) {  // check if Y scale is nonlinear, and if so, do transform
+		 if(!lowerSlopeNatural) lowerSlope *= YY.fInPrime(func(xmin));
+		 if(!upperSlopeNatural) upperSlope *= YY.fInPrime(func(xmax));
+	 }
+	 // note that each of ends has 3 points with two equal gaps, since they were obtained by bisection
+	 // so the step sizes are easy to get
+	 // the 'natural slope' option for sampled functions has a different meaning than 
+	 // for normal splines.  In this case, the derivative is adjusted to make the
+	 // second derivative constant on the last two points at each end
+	 // which is consistent with the error sampling technique we used to get here
+	 if(lowerSlopeNatural) {
+		 float_type hlower=X[1]-X[0];
+		 lowerSlope=0.5*(-F[2]-3*F[0]+4*F[1])/hlower;
+		 lowerSlopeNatural=false; // it's not the usual meaning of natural any more 
+	 }
+	 if(upperSlopeNatural) {
+		 float_type hupper=X[np-1]-X[np-2];
+		 upperSlope=0.5*(F[np-3]+3*F[np-1]-4*F[np-2])/hupper;
+		 upperSlopeNatural=false; // it's not the usual meaning of natural any more 
+	 }
+	 this->set_domain(xmin, xmax);
+	 
+	 spline(lowerSlopeNatural, lowerSlope, upperSlopeNatural, upperSlope);
+	 lastKLow=0;
+	 keepme.release_for_return();
+	 return *this;
 }
 
 //  This function is the reason for this class to exist
 // it computes the interpolated function, and (if requested) its proper first and second derivatives including all coordinate transforms
-template <typename float_type> float_type interpolating_function<float_type>::value_with_derivatives(
+template <typename float_type> float_type interpolating_function_p<float_type>::value_with_derivatives(
 				float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
 {
+	if(sampler_function.valid()) { 
+		// if this is non-null, we are sampling data for later, so just return raw function
+		// however, transform it into our sampling space, first.
+		if(yprime) *yprime=0;
+		if(yprime2) *yprime2=0;
+		sampler_function->increment_evaluations();
+		return fTransform.Y.fIn(sampler_function(fTransform.X.fOut(x))); // derivatives are completely undefined
+	}
+					
 	if(x < this->xmin() || x > this->xmax()) {
 		std::ostringstream outstr;
@@ -613,15 +918,19 @@
 	float_type xraw=x;
 	
-	// template here is impossible! if(fXin && fXin != (Identity<float_type>) ) 
-	x=fXin(x); // save time by explicitly testing for identity function here
+	if(fTransform.X.fTransformed) x=fTransform.X.fHasStaticTransforms? 
+		fTransform.X.pIn(x) : fTransform.X.fIn(x); // save time by explicitly testing for identity function here
 	
 	int klo=0, khi=X.size()-1;
+
+	if(khi < 0) throw c2_exception("Uninitialized interpolating function being evaluated");	
+	
+	const float_type *XX=&X[lastKLow]; // make all fast checks short offsets from here
 	
 	if(!xInverted) { // select search depending on whether transformed X is increasing or decreasing
-		if(lastKLow >=0 && (X[lastKLow] <= x) && (X[lastKLow+1] >= x) ) { // already bracketed
+		if((XX[0] <= x) && (XX[1] >= x) ) { // already bracketed
 			klo=lastKLow;
-		} else if(lastKLow >=0 && (X[lastKLow+1] <= x) && (X[lastKLow+2] > x)) { // in next bracket to the right
+		} else if((XX[1] <= x) && (XX[2] >= x)) { // in next bracket to the right
 			klo=lastKLow+1;		
-		} else if(lastKLow > 0 && (X[lastKLow-1] <= x) && (X[lastKLow] > x)) { // in next bracket to the left
+		} else if(lastKLow > 0 && (XX[-1] <= x) && (XX[0] >= x)) { // in next bracket to the left
 			klo=lastKLow-1;		
 		} else { // not bracketed, not close, start over
@@ -634,9 +943,9 @@
 		}
 	} else {
-		if(lastKLow >=0 && (X[lastKLow] >= x) && (X[lastKLow+1] <= x) ) { // already bracketed
+		if((XX[0] >= x) && (XX[1] <= x) ) { // already bracketed
 			klo=lastKLow;
-		} else if(lastKLow >=0 && (X[lastKLow+1] >= x) && (X[lastKLow+2] < x)) { // in next bracket to the right
+		} else if((XX[1] >= x) && (XX[2] <= x)) { // in next bracket to the right
 			klo=lastKLow+1;		
-		} else if(lastKLow > 0 && (X[lastKLow-1] >= x) && (X[lastKLow] < x)) { // in next bracket to the left
+		} else if(lastKLow > 0 && (XX[-1] >= x) && (XX[0] <= x)) { // in next bracket to the left
 			klo=lastKLow-1;		
 		} else { // not bracketed, not close, start over
@@ -659,32 +968,25 @@
 	float_type ylo=F[klo], yhi=F[khi], y2lo=y2[klo], y2hi=y2[khi];
 	float_type y=a*ylo+b*yhi+((a*a*a-a)*y2lo+(b*b*b-b)*y2hi)*(h*h)/6.0;
-	
-	// template here is impossible! if(fYin && fYin != Identity) 
-	y=fYout(y); // save time by explicitly testing for identity function here
+
+	float_type yp0=0; // the derivative in interpolating table coordinates
+	float_type ypp0=0; // second derivative
 	
 	if(yprime || yprime2) {
-		float_type fpi=1.0/fYinPrime(y);
-		float_type gp=fXinPrime(xraw);
-		float_type yp0=(yhi-ylo)/h+((3*b*b-1)*y2hi-(3*a*a-1)*y2lo)*h/6.0; // the derivative in interpolating table coordinates
-		
-		// from Mathematica Dt[InverseFunction[f][y[g[x]]], x]
-		if(yprime) *yprime=gp*yp0*fpi; // the real derivative of the inverse transformed output 
-		if(yprime2) {
-			float_type ypp0=b*y2hi+a*y2lo;
-			float_type fpp=fYinDPrime(y);
-			float_type gpp=fXinDPrime(xraw);
-			// also from Mathematica Dt[InverseFunction[f][y[g[x]]], {x,2}]
-			if(yprime2) *yprime2=(gp*gp*ypp0 + yp0*gpp - gp*gp*yp0*yp0*fpp*fpi*fpi)*fpi; 
-		}
-	}
-	
-	return y;
-}
-
-template <typename float_type> void interpolating_function<float_type>::set_lower_extrapolation(float_type bound) 
+		yp0=(yhi-ylo)/h+((3*b*b-1)*y2hi-(3*a*a-1)*y2lo)*h/6.0; // the derivative in interpolating table coordinates
+		ypp0=b*y2hi+a*y2lo; // second derivative
+	}
+	
+	if(fTransform.isIdentity) {
+		if(yprime) *yprime=yp0;
+		if(yprime2) *yprime2=ypp0;
+		return y;	
+	} else return fTransform.evaluate(xraw, y, yp0, ypp0, yprime, yprime2);
+}
+
+template <typename float_type> void interpolating_function_p<float_type>::set_lower_extrapolation(float_type bound) 
 {
 	int kl = 0 ;
 	int kh=kl+1;
-	float_type xx=fXin(bound);
+	float_type xx=fTransform.X.fIn(bound);
 	float_type h0=X[kh]-X[kl];
 	float_type h1=xx-X[kl];
@@ -702,9 +1004,9 @@
 }
 
-template <typename float_type> void interpolating_function<float_type>::set_upper_extrapolation(float_type bound) 
+template <typename float_type> void interpolating_function_p<float_type>::set_upper_extrapolation(float_type bound) 
 {
 	int kl = X.size()-2 ;
 	int kh=kl+1;
-	float_type xx=fXin(bound);
+	float_type xx=fTransform.X.fIn(bound);
 	float_type h0=X[kh]-X[kl];
 	float_type h1=xx-X[kl];
@@ -721,18 +1023,4 @@
 }
 
-// move derivatives into our internal coordinates (use splint to go the other way!)
-template <typename float_type> void interpolating_function<float_type>::localize_derivatives(
-   float_type xraw, float_type y, float_type yp, float_type ypp, float_type *y0, float_type *yprime, float_type *yprime2) const 
-{
-	float_type fp=fYinPrime(y);
-	float_type gp=fXinPrime(xraw);
-	float_type fpp=fYinDPrime(y);
-	float_type gpp=fXinDPrime(xraw);
-	
-	if(y0) *y0=fYin(y);
-	if(yprime) *yprime=yp*fp/gp; // Mathematica Dt[f[y[InverseFunction[g][x]]], x]
-	if(yprime2) *yprime2=( yp*yp*fpp - fp*yp*gpp/gp + fp*ypp )/(gp*gp) ; // Mathematica Dt[f[y[InverseFunction[g][x]]], {x,2}]
-}
-
 // return a new interpolating_function which is the unary function of an existing interpolating_function
 // can also be used to generate a resampling of another c2_function on a different grid
@@ -740,26 +1028,27 @@
 // and doing b=a.unary_operator(c) where c is a c2_function (probably another interpolating_function)
 
-template <typename float_type> interpolating_function<float_type>& 
-	interpolating_function<float_type>::unary_operator(const c2_function<float_type> &source) const 
+template <typename float_type> interpolating_function_p<float_type>& 
+	interpolating_function_p<float_type>::unary_operator(const c2_function<float_type> &source) const 
 {
 	size_t np=X.size();
 	std::vector<float_type>yv(np);
-	c2_composed_function<float_type> comp(source, *this);
+	c2_ptr<float_type> comp(source(*this));
 	float_type yp0, yp1, ypp;
 	
-	for(size_t i=0; i<np; i++) { 
-		yv[i]=source(fYout(F[i])); // copy pointwise the function of our data values
-	}
-		
-	comp(Xraw.front(), &yp0, &ypp); // get derivative at front
-	comp(Xraw.back(), &yp1, &ypp); // get derivative at back
-	
-	return  *new interpolating_function(Xraw, yv, false, yp0, false, yp1,
-									   fXin, fXinPrime, fXinDPrime,
-									   fYin, fYinPrime, fYinDPrime, fYout);
+	for(size_t i=1; i<np-1; i++) { 
+		yv[i]=source(fTransform.Y.fOut(F[i])); // copy pointwise the function of our data values
+	}
+		
+	yv.front()=comp(Xraw.front(), &yp0, &ypp); // get derivative at front
+	yv.back()= comp(Xraw.back(), &yp1, &ypp); // get derivative at back
+	
+	interpolating_function_p &copy=clone();
+	copy.load(this->Xraw, yv, false, yp0, false, yp1);
+	
+	return copy;
 }
 
 template <typename float_type> void 
-interpolating_function<float_type>::get_data(std::vector<float_type> &xvals, std::vector<float_type> &yvals) const throw()
+interpolating_function_p<float_type>::get_data(std::vector<float_type> &xvals, std::vector<float_type> &yvals) const throw()
 {
 	
@@ -767,88 +1056,34 @@
 	yvals.resize(F.size());
 	
-	for(size_t i=0; i<F.size(); i++) yvals[i]=fYout(F[i]);
-}
-
-template <typename float_type> interpolating_function<float_type> & 
-	interpolating_function<float_type>::binary_operator(const c2_function<float_type> &rhs,
-		c2_binary_function<float_type> *combining_stub) const
+	for(size_t i=0; i<F.size(); i++) yvals[i]=fTransform.Y.fOut(F[i]);
+}
+
+template <typename float_type> interpolating_function_p<float_type> & 
+	interpolating_function_p<float_type>::binary_operator(const c2_function<float_type> &rhs,
+		const c2_binary_function<float_type> *combining_stub) const
 {	
 	size_t np=X.size();	
 	std::vector<float_type> yv(np);
-	c2_constant<float_type> fval;
-	c2_constant<float_type> yval;
+	c2_constant_p<float_type> fval(0);
 	float_type yp0, yp1, ypp;
 	
-	for(size_t i=0; i<np; i++) { 
-		fval.reset(fYout(F[i])); // update the constant function pointwise
-		yval.reset(rhs(Xraw[i])); 
-		yv[i]=(*combining_stub).combine(fval, yval, Xraw[i], (float_type *)0, (float_type *)0); // compute rhs & combine without derivatives
-	}
-	
-	(*combining_stub).combine(*this, rhs, Xraw.front(), &yp0, &ypp); // get derivative at front
-	(*combining_stub).combine(*this, rhs, Xraw.back(),  &yp1, &ypp); // get derivative at back
-
-	delete combining_stub;
-	
-	return  *new interpolating_function(Xraw, yv, false, yp0, false, yp1,
-									   fXin, fXinPrime, fXinDPrime,
-									   fYin, fYinPrime, fYinDPrime, fYout);	
-}
-
-template <typename float_type>  float_type c2_f_logprime(float_type x) { return 1.0/x; } // the derivative of log(x)
-template <typename float_type>  float_type c2_f_logprime2(float_type x) { return -1.0/(x*x); } // the second derivative of log(x)
-
-template <typename float_type> log_lin_interpolating_function<float_type>::log_lin_interpolating_function(
-										const std::vector<float_type> &x, const std::vector<float_type> &f, 
-										bool lowerSlopeNatural, float_type lowerSlope, 
-										bool upperSlopeNatural, float_type upperSlope) 
-										: interpolating_function<float_type>()
-{
-	init(x, f, lowerSlopeNatural, lowerSlope, upperSlopeNatural, upperSlope,
-		 (float_type (*)(float_type)) (std::log), c2_f_logprime, c2_f_logprime2, 0, 0, 0, 0);
-}
-
-template <typename float_type> lin_log_interpolating_function<float_type>::lin_log_interpolating_function(
-										const std::vector<float_type> &x, const std::vector<float_type> &f, 
-										bool lowerSlopeNatural, float_type lowerSlope, 
-										bool upperSlopeNatural, float_type upperSlope)
-										: interpolating_function<float_type>()
-{
-	init(x, f, lowerSlopeNatural, lowerSlope, upperSlopeNatural, upperSlope,
-		 0, 0, 0, 
-		 (float_type (*)(float_type)) (std::log), c2_f_logprime, c2_f_logprime2, 
-		 (float_type (*)(float_type)) (std::exp) );
-}
-
-template <typename float_type> log_log_interpolating_function<float_type>::log_log_interpolating_function(
-										const std::vector<float_type> &x, const std::vector<float_type> &f, 
-										bool lowerSlopeNatural, float_type lowerSlope, 
-										bool upperSlopeNatural, float_type upperSlope) 
-										: interpolating_function<float_type>()
-{
-	init(x, f, lowerSlopeNatural, lowerSlope, upperSlopeNatural, upperSlope,
-		 (float_type (*)(float_type)) (std::log), c2_f_logprime, c2_f_logprime2, 
-		 (float_type (*)(float_type)) (std::log), c2_f_logprime, c2_f_logprime2, 
-		 (float_type (*)(float_type)) (std::exp) );
-}
-
-template <typename float_type> float_type c2_f_recip(float_type x) { return 1.0/x; }
-template <typename float_type> float_type c2_f_recipprime(float_type x) { return -1.0/(x*x); } // the derivative of 1/x
-template <typename float_type> float_type c2_f_recipprime2(float_type x) { return 2.0/(x*x*x); } // the second derivative of 1/x
-
-template <typename float_type> arrhenius_interpolating_function<float_type>::arrhenius_interpolating_function(
-										const std::vector<float_type> &x, const std::vector<float_type> &f, 
-										bool lowerSlopeNatural, float_type lowerSlope, 
-										bool upperSlopeNatural, float_type upperSlope) 
-										: interpolating_function<float_type>()
-{
-	init(x, f, lowerSlopeNatural, lowerSlope, upperSlopeNatural, upperSlope,
-		c2_f_recip, c2_f_recipprime, c2_f_recipprime2, 
-		(float_type (*)(float_type)) (std::log), c2_f_logprime, c2_f_logprime2, 
-		(float_type (*)(float_type)) (std::exp) );
-}
-
-template <typename float_type> c2_inverse_function<float_type>::c2_inverse_function(const c2_function<float_type> &source) 
-	: c2_plugin_function<float_type>(source) 
+	c2_const_ptr<float_type> stub(*combining_stub);  // manage ownership
+	
+	for(size_t i=1; i<np-1; i++) { 
+		fval.reset(fTransform.Y.fOut(F[i])); // update the constant function pointwise
+		yv[i]=combining_stub->combine(fval, rhs, Xraw[i], (float_type *)0, (float_type *)0); // compute rhs & combine without derivatives
+	}
+	
+	yv.front()=combining_stub->combine(*this, rhs, Xraw.front(), &yp0, &ypp); // get derivative at front
+	yv.back()= combining_stub->combine(*this, rhs, Xraw.back(),  &yp1, &ypp); // get derivative at back
+	
+	interpolating_function_p &copy=clone();
+	copy.load(this->Xraw, yv, false, yp0, false, yp1);
+	
+	return copy;
+}
+
+template <typename float_type> c2_inverse_function_p<float_type>::c2_inverse_function_p(const c2_function<float_type> &source) 
+	: c2_function<float_type>(), func(source)
 {
 	float_type l=source.xmin();
@@ -864,5 +1099,5 @@
 }
 
-template <typename float_type> float_type c2_inverse_function<float_type>::value_with_derivatives(
+template <typename float_type> float_type c2_inverse_function_p<float_type>::value_with_derivatives(
 					float_type x, float_type *yprime, float_type *yprime2
 			) const throw(c2_exception)
@@ -927,10 +1162,10 @@
 		for(int i=0; i<=np; i++) cum[i]*=m;
 	}
-	if(inverse_function) interpolating_function<float_type>(cum, be); // use cum as x axis in inverse function
-	else interpolating_function<float_type>(be, cum); // else use lower bin edge as x axis
+	if(inverse_function) interpolating_function_p<float_type>(cum, be); // use cum as x axis in inverse function
+	else interpolating_function_p<float_type>(be, cum); // else use lower bin edge as x axis
 	std::fill(this->y2.begin(), this->y2.end(), 0.0); // clear second derivatives, to we are piecewise linear
 }
 
-template <typename float_type> c2_piecewise_function<float_type>::c2_piecewise_function() 
+template <typename float_type> c2_piecewise_function_p<float_type>::c2_piecewise_function_p() 
 : c2_function<float_type>(), lastKLow(-1)
 {
@@ -938,11 +1173,9 @@
 }
 
-template <typename float_type> c2_piecewise_function<float_type>::~c2_piecewise_function() 
-{
-	size_t np=functions.size();
-	for(size_t i=0; i<np; i++) if(owns[i]) delete functions[i];
-}
-
-template <typename float_type> float_type c2_piecewise_function<float_type>::value_with_derivatives(
+template <typename float_type> c2_piecewise_function_p<float_type>::~c2_piecewise_function_p() 
+{
+}
+
+template <typename float_type> float_type c2_piecewise_function_p<float_type>::value_with_derivatives(
 		  float_type x, float_type *yprime, float_type *yprime2
 		  ) const throw(c2_exception)
@@ -974,9 +1207,10 @@
 }
 
-template <typename float_type> void c2_piecewise_function<float_type>::append_function(
-	c2_function<float_type> &func, bool pass_ownership) throw(c2_exception)
-{
+template <typename float_type> void c2_piecewise_function_p<float_type>::append_function(
+	const c2_function<float_type> &func) throw(c2_exception)
+{
+	c2_const_ptr<float_type> keepfunc(func);  //  manage function before we can throw any exceptions
 	if(functions.size()) { // check whether there are any gaps to fill, etc.
-		c2_function<float_type> &tail=*(functions.back());
+		const c2_function<float_type> &tail=functions.back();
 		float_type x0=tail.xmax();
 		float_type x1=func.xmin();
@@ -985,78 +1219,333 @@
 			float_type y0=tail(x0);
 			float_type y1=func(x1);
-			c2_function<float_type> *connector=new c2_linear<float_type>(x0, y0, (y1-y0)/(x1-x0));
-			connector->set_domain(x0,x1);
-			functions.push_back(connector);
-			owns.push_back(true);
+			c2_function<float_type> &connector=*new c2_linear_p<float_type>(x0, y0, (y1-y0)/(x1-x0));
+			connector.set_domain(x0,x1);
+			functions.push_back(c2_const_ptr<float_type>(connector));
 			this->sampling_grid->push_back(x1);
 		} else if(x0>x1) throw c2_exception("function domains not increasing in c2_piecewise_function");
 	}
-	functions.push_back(&func);
-	owns.push_back(pass_ownership);
+	functions.push_back(keepfunc);
 	// extend our domain to include all known functions
 	this->set_domain(functions.front()->xmin(), functions.back()->xmax());
 	// extend our sampling grid with the new function's grid, with the first point dropped to avoid duplicates
-	std::vector<float_type> &newgrid=func.get_sampling_grid(func.xmin(), func.xmax());
+	std::vector<float_type> newgrid;
+	func.get_sampling_grid(func.xmin(), func.xmax(), newgrid);
 	this->sampling_grid->insert(this->sampling_grid->end(), newgrid.begin()+1, newgrid.end());
-	delete &newgrid;
-}
-
-template <typename float_type> c2_connector_function<float_type>::c2_connector_function(
-	const c2_function<float_type> &f1, const c2_function<float_type> &f2, float_type x0, float_type x2, 
+}
+
+template <typename float_type> c2_connector_function_p<float_type>::c2_connector_function_p(
+	float_type x0, const c2_function<float_type> &f0, float_type x2, const c2_function<float_type> &f2, 
 	bool auto_center, float_type y1)
-
-: c2_function<float_type>()
-{
-	float_type y0, yp0, ypp0, y2, yp2, ypp2;
-	fdx=(x2-x0)/2.0;
-	fhinv=1.0/fdx;
-	fx1=(x0+x2)/2.0;
-	
-	y0=f1.value_with_derivatives(x0, &yp0, &ypp0); // get left wall values from conventional computation
-	y2=f2.value_with_derivatives(x2, &yp2, &ypp2); // get right wall values from conventional computation
+	: c2_function<float_type>()
+{
+	c2_const_ptr<float_type> left(f0), right(f2); // make sure if these are unowned, they get deleted
+	c2_fblock<float_type> fb0, fb2;
+	fb0.x=x0;
+	f0.fill_fblock(fb0);
+	fb2.x=x2;
+	f2.fill_fblock(fb2);
+	init(fb0, fb2, auto_center, y1);
+}
+
+template <typename float_type> c2_connector_function_p<float_type>::c2_connector_function_p(
+	float_type x0, float_type y0, float_type yp0, float_type ypp0, 
+	float_type x2, float_type y2, float_type yp2, float_type ypp2, 
+	bool auto_center, float_type y1)
+	: c2_function<float_type>()
+{
+	c2_fblock<float_type> fb0, fb2;
+	fb0.x=x0; fb0.y=y0; fb0.yp=yp0; fb0.ypp=ypp0;
+	fb2.x=x2; fb2.y=y2; fb2.yp=yp2; fb2.ypp=ypp2;	
+	init(fb0, fb2, auto_center, y1);
+}
+
+template <typename float_type> c2_connector_function_p<float_type>::c2_connector_function_p(
+	const c2_fblock<float_type> &fb0,
+	const c2_fblock<float_type> &fb2,
+	bool auto_center, float_type y1)
+	: c2_function<float_type>()
+{
+	init(fb0, fb2, auto_center, y1);
+}
+
+template <typename float_type> void c2_connector_function_p<float_type>::init(
+	const c2_fblock<float_type> &fb0,
+	const c2_fblock<float_type> &fb2,
+	bool auto_center, float_type y1)
+{
+	float_type dx=(fb2.x-fb0.x)/2.0;
+	fhinv=1.0/dx;
 
 	// scale derivs to put function on [-1,1] since mma  solution is done this way
-	yp0*=fdx;
-	yp2*=fdx;
-	ypp0*=fdx*fdx;
-	ypp2*=fdx*fdx;
-	
-	float_type ff0=(8*(y0 + y2) + 5*(yp0 - yp2) + ypp0 + ypp2)/16.0;
+	float_type yp0=fb0.yp*dx;
+	float_type yp2=fb2.yp*dx;
+	float_type ypp0=fb0.ypp*dx*dx;
+	float_type ypp2=fb2.ypp*dx*dx;
+	
+	float_type ff0=(8*(fb0.y + fb2.y) + 5*(yp0 - yp2) + ypp0 + ypp2)*0.0625;
 	if(auto_center) y1=ff0; // forces ff to be 0 if we are auto-centering
 	
-	// y[x_] = y1 + x (a + b x) + (x-1) x (x+1) (c + d x + e x^2 + f x^3)
+	// y[x_] = y1 + x (a + b x) + x [(x-1) (x+1)] (c + d x) + x (x-1)^2 (x+1)^2 (e + f x)
+	// y' = a + 2 b x + d x [(x+1)(x-1)] + (c + d x)(3x^2-1) + f x [(x+1)(x-1)]^2 + (e + f x)[(x+1)(x-1)](5x^2-1)  
+	// y'' = 2 b + 6x(c + d x) + 2d(3x^2-1) + 4x(e + f x)(5x^2-3) + 2f(x^2-1)(5x^2-1)
 	fy1=y1; 
-	fa=-(y0 - y2)/2.;
-	fb=(y0 - 2*y1 + y2)/2.;
-	fc=(7*(y0 - y2 + yp0 + yp2) + ypp0 - ypp2)/16.;
-	fd=(32*y1 - 16*(y2 + y0) + 9*(yp2 - yp0) - ypp0 - ypp2)/16.;
-	fe=(3*(y2 - y0 - yp0 - yp2) - ypp0 + ypp2)/16.;
+	fa=(fb2.y - fb0.y)*0.5;
+	fb=(fb0.y + fb2.y)*0.5 - y1;
+	fc=(yp2+yp0-2.*fa)*0.25;
+	fd=(yp2-yp0-4.*fb)*0.25;
+	fe=(ypp2-ypp0-12.*fc)*0.0625;
 	ff=(ff0 - y1);
-	// y'[x] = a + 2 b x + (3x^2 - 1)   (c + d x + e x^2 + f x^3) + (x-1) x (x+1) (d + 2 e x + 3 f x^2 )
-	// y''[x] = 2b + (x-1) x (x+1) (2 e + 6 f x) + 2 (3 x^2 -1) (d + 2 e x + 3 f x^2 ) + 6 x (c + d x + e x^2 + f x^3)
-	this->set_domain(x0,x2); // this is where the function is valid
-}
-
-template <typename float_type> c2_connector_function<float_type>::~c2_connector_function() 
-{
-}
-
-template <typename float_type> float_type c2_connector_function<float_type>::value_with_derivatives(
+	this->set_domain(fb0.x, fb2.x); // this is where the function is valid
+}
+
+template <typename float_type> c2_connector_function_p<float_type>::~c2_connector_function_p() 
+{
+}
+
+template <typename float_type> float_type c2_connector_function_p<float_type>::value_with_derivatives(
 	float_type x, float_type *yprime, float_type *yprime2
 	) const throw(c2_exception)
 {
-																									
-	float_type dx=(x-fx1)*fhinv;
-	float_type q1=fc + dx*(fd + dx*(fe + dx*ff));
-	float_type xp1=(dx-1)*(dx+1)*dx;
-	
-	float_type y= fy1 + dx*(fa+fb*dx) + xp1*q1;
+	float_type x0=this->xmin(), x2=this->xmax();
+	float_type dx=(x-(x0+x2)*0.5)*fhinv;
+	float_type q1=(x-x0)*(x-x2)*fhinv*fhinv; // exactly vanish all bits at both ends
+	float_type q2=dx*q1;
+	
+	float_type r1=fa+fb*dx;
+	float_type r2=fc+fd*dx;
+	float_type r3=fe+ff*dx;
+	
+	float_type y=fy1+dx*r1+q2*r2+q1*q2*r3;
+	
 	if(yprime || yprime2) {
-		float_type q2 =fd + dx*(2*fe + dx*3*ff);	
-		float_type q3=2*fe+6*ff*dx;
-		float_type xp2=(3*dx*dx-1);
-		if(yprime) *yprime=(fa + 2*fb*dx + xp2*q1 + xp1*q2)*fhinv;
-		if(yprime2) *yprime2=(2*fb+xp1*q3+2*xp2*q2+6*dx*q1)*fhinv*fhinv;
+		float_type q3=3*q1+2;
+		float_type q4=5*q1+4;
+		if(yprime) *yprime=(fa+2*fb*dx+fd*q2+r2*q3+ff*q1*q2+q1*q4*r3)*fhinv;
+		if(yprime2) *yprime2=2*(fb+fd*q3+3*dx*r2+ff*q1*q4+r3*(2*dx*(5*q1+2)))*fhinv*fhinv;
 	}
 	return y;
 }
+
+// the recursive part of the sampler is agressively designed to minimize copying of data... lots of pointers
+template <typename float_type> void c2_function<float_type>::sample_step(c2_sample_recur &rb) const throw(c2_exception)
+{
+	std::vector< recur_item > &rb_stack=*rb.rb_stack; // heap-based stack of data for recursion
+	rb_stack.clear();
+	
+	recur_item top;
+	top.depth=0; top.done=false; top.f0index=0; top.f2index=0;
+	
+	// push storage for our initial elements
+	rb_stack.push_back(top);
+	rb_stack.back().f1=*rb.f0;
+	rb_stack.back().done=true;
+
+	rb_stack.push_back(top);
+	rb_stack.back().f1=*rb.f1;
+	rb_stack.back().done=true;
+	
+	if(!rb.inited) {
+		rb.dx_tolerance=10.0*std::numeric_limits<float_type>::epsilon();
+		rb.abs_tol_min=10.0*std::numeric_limits<float_type>::min();
+		rb.inited=true;
+	}
+
+	// now, push our first real element
+	top.f0index=0; // left element is stack[0]
+	top.f2index=1; // right element is stack[1]
+	rb_stack.push_back(top);
+
+	while(rb_stack.size() > 2) {
+		recur_item &back=rb_stack.back();
+		if(back.done) {
+			rb_stack.pop_back();
+			continue;
+		}
+		back.done=true;
+
+		c2_fblock<float_type> &f0=rb_stack[back.f0index].f1, &f2=rb_stack[back.f2index].f1; 
+		c2_fblock<float_type> &f1=back.f1; // will hold new middle values
+		size_t f1index=rb_stack.size()-1; // our current offset
+
+		// std::cout << "processing: " << rb_stack.size() << " " << 
+		// (&back-&rb_stack.front())  << " " << back.depth << " " << f0.x << " " << f2.x << std::endl;
+
+		f1.x=0.5*(f0.x + f2.x); // center of interval
+		float_type dx2=0.5*(f2.x - f0.x);
+			
+		// check for underflow on step size, which prevents us from achieving specified accuracy. 
+		if(std::abs(dx2) < std::abs(f1.x)*rb.dx_tolerance || std::abs(dx2) < rb.abs_tol_min) {
+			std::ostringstream outstr;
+			outstr << "Step size underflow in adaptive_sampling at depth=" << back.depth << ", x= " << f1.x;
+			throw c2_exception(outstr.str().c_str());
+		}
+		
+		fill_fblock(f1);
+		
+		if(c2_isnan(f1.y) || f1.ypbad || f1.yppbad) {
+			// can't go any further if a nan has appeared
+			bad_x_point=f1.x;
+			throw c2_exception("NaN encountered while sampling function"); 
+		}
+
+		float_type eps; 
+		if(rb.derivs==2) {
+			// this is code from connector_function to compute the value at the midpoint
+			// it is re-included here to avoid constructing a complete c2connector
+			// just to find out if we are close enough
+			float_type ff0=(8*(f0.y + f2.y) + 5*(f0.yp - f2.yp)*dx2 + (f0.ypp+f2.ypp)*dx2*dx2)*0.0625;
+			// we are converging as at least x**5 and bisecting, so real error on final step is smaller
+			eps=std::abs(ff0-f1.y)/32.0;   
+		} else {
+			// there are two tolerances to meet... the shift in the estimate of the actual point,
+			// and the difference between the current points and the extremum
+			// build all the coefficients needed to construct the local parabola
+			float_type ypcenter, ypp;
+			if (rb.derivs==1) {
+				// linear extrapolation error using exact derivs
+				eps = (std::abs(f0.y+f0.yp*dx2-f1.y)+std::abs(f2.y-f2.yp*dx2-f1.y))*0.125; 
+				ypcenter=2*f1.yp*dx2; // first deriv scaled so this interval is on [-1,1]
+				ypp=2*(f2.yp-f0.yp)*dx2*dx2; // second deriv estimate scaled so this interval is on [-1,1]
+			} else {
+				// linear interpolation error without derivs if we are at top level
+				//  or 3-point parabolic interpolation estimates from previous level, if available
+				ypcenter=(f2.y-f0.y)*0.5; // derivative estimate at center
+				ypp=(f2.y+f0.y-2*f1.y); // second deriv estimate
+				if(back.depth==0) eps=std::abs((f0.y+f2.y)*0.5 - f1.y)*2; // penalize first step
+				else eps=std::abs(f1.y-back.previous_estimate)*0.25; 
+			}
+			float_type ypleft=ypcenter-ypp; // derivative at left edge
+			float_type ypright=ypcenter+ypp; // derivative at right edge
+			float_type extremum_eps=0;
+			if((ypleft*ypright) <=0) // y' changes sign if we have an extremum
+			{
+				// compute position and value of the extremum this way
+				float_type xext=-ypcenter/ypp;
+				float_type yext=f1.y + xext*ypcenter + 0.5*xext*xext*ypp;
+				// and then find the the smallest offset of it from a point, looking in the left or right side
+				if(xext <=0) extremum_eps=std::min(std::abs(f0.y-yext), std::abs(f1.y-yext));
+				else extremum_eps=std::min(std::abs(f2.y-yext), std::abs(f1.y-yext));
+			}
+			eps=std::max(eps, extremum_eps); // if previous shot was really bad, keep trying
+		}
+		
+		if(eps < rb.abs_tol ||  eps < std::abs(f1.y)*rb.rel_tol) {
+			if(rb.out) {
+				// we've met the tolerance, and are building a function, append two connectors
+				rb.out->append_function(
+					*new c2_connector_function_p<float_type>(f0, f1, true, 0.0)
+				);
+				rb.out->append_function(
+					*new c2_connector_function_p<float_type>(f1, f2, true, 0.0)
+				);
+			}
+			if(rb.xvals && rb.yvals) {
+				rb.xvals->push_back(f0.x);
+				rb.xvals->push_back(f1.x);
+				rb.yvals->push_back(f0.y);
+				rb.yvals->push_back(f1.y);
+				// the value at f2 will get pushed in the next segment... it is not forgotten
+			}
+		} else {
+			top.depth=back.depth+1; // increment depth counter
+
+			// save the last things we need from back before a push happens, in case 
+			// the push causes a reallocation and moves the whole stack.
+			size_t f0index=back.f0index, f2index=back.f2index;
+			float_type left=0, right=0;
+			if(rb.derivs==0) {
+				// compute three-point parabolic interpolation estimate of right-hand and left-hand midpoint
+				left=(6*f1.y + 3*f0.y - f2.y) * 0.125;
+				right=(6*f1.y + 3*f2.y - f0.y) * 0.125; 
+			}
+			
+			top.f0index=f1index; top.f2index=f2index; // insert pointers to right side data into our recursion block			
+			top.previous_estimate=right; 
+			rb_stack.push_back(top);
+
+			top.f0index=f0index; top.f2index=f1index; // insert pointers to left side data into our recursion block
+			top.previous_estimate=left; 
+			rb_stack.push_back(top);
+		}
+	}
+}
+
+template <typename float_type>  c2_piecewise_function_p<float_type> *
+	c2_function<float_type>::adaptively_sample(
+		float_type xmin, float_type xmax,
+		float_type abs_tol, float_type rel_tol,
+		int derivs, std::vector<float_type> *xvals, std::vector<float_type> *yvals) const throw(c2_exception)
+{
+	c2_fblock<float_type> f0, f2;
+	c2_sample_recur rb;
+	std::vector< recur_item > rb_stack;
+	rb_stack.reserve(20); // enough for most operations
+	rb.rb_stack=&rb_stack;
+	rb.out=0;
+	if(derivs==2) rb.out=new c2_piecewise_function_p<float_type>();
+	c2_ptr<float_type> pieces(*rb.out); // manage this function, if any, so it deletes on an exception
+	rb.rel_tol=rel_tol;
+	rb.abs_tol=abs_tol;
+	rb.xvals=xvals;
+	rb.yvals=yvals;
+	rb.derivs=derivs;
+	rb.inited=false;
+	
+	if(xvals && yvals) {
+		xvals->clear();
+		yvals->clear();
+	}
+	 
+	// create xgrid as a automatic-variable copy of the sampling grid so the exception handler correctly 
+	// disposes of it.
+	std::vector<float_type> xgrid;
+	get_sampling_grid(xmin, xmax, xgrid);
+	int np=xgrid.size();	 
+
+	f2.x=xgrid[0];
+	fill_fblock(f2);
+	if(c2_isnan(f2.y) || f2.ypbad || f2.yppbad) {
+		// can't go any further if a nan has appeared
+		bad_x_point=f2.x;
+		throw c2_exception("NaN encountered while sampling function"); 
+	}
+
+	for(int i=0; i<np-1; i++) {
+		f0=f2; // copy upper bound to lower before computing new upper bound
+
+		f2.x=xgrid[i+1];
+		fill_fblock(f2);
+		if(c2_isnan(f2.y) || f2.ypbad || f2.yppbad) {
+			// can't go any further if a nan has appeared
+			bad_x_point=f2.x;
+			throw c2_exception("NaN encountered while sampling function"); 
+		}
+
+		rb.f0=&f0; rb.f1=&f2;
+		sample_step(rb);
+	}
+	if(xvals && yvals) { // push final point in vector
+		xvals->push_back(f2.x);
+		yvals->push_back(f2.y);
+	}
+	
+	if(rb.out) rb.out->set_sampling_grid(xgrid); // reflect old sampling grid, which still should be right
+	pieces.release_for_return(); // unmanage the piecewise_function so we can return it
+	return rb.out;
+}
+
+template <typename float_type, typename Final> 
+ interpolating_function_p<float_type> & inverse_integrated_density_function(
+	const std::vector<float_type> &bincenters, const c2_function<float_type> &binheights)
+	throw(c2_exception)
+{	
+	return (new Final())->load_random_generator_function(bincenters, binheights);
+}
+
+template <typename float_type, typename Final> 
+	interpolating_function_p<float_type> & inverse_integrated_density_bins(
+	const std::vector<float_type> &bins, const std::vector<float_type> &binheights)
+	throw(c2_exception)
+{	
+	return (new Final())->load_random_generator_bins(bins, binheights);
+}
