Changeset 1230 for trunk/examples/extended/electromagnetic/TestEm7/include

Timestamp:

Jan 8, 2010, 3:02:48 PM (14 years ago)

Author:

garnier

Message:

update to geant4.9.3

Location:

trunk/examples/extended/electromagnetic/TestEm7/include

Files:

• DetectorConstruction.hh (modified) (1 diff)
• DetectorMessenger.hh (modified) (1 diff)
• EventAction.hh (modified) (1 diff)
• EventActionMessenger.hh (modified) (1 diff)
• G4ScreenedNuclearRecoil.hh (modified) (16 diffs)
• PhysListEmLivermore.hh (modified) (1 diff)
• PhysListEmPenelope.hh (modified) (1 diff)
• PhysListEmStandard.hh (modified) (1 diff)
• PhysListEmStandardNR.hh (modified) (1 diff)
• PhysListEmStandardSS.hh (modified) (1 diff)
• PhysicsList.hh (modified) (5 diffs)
• PhysicsListMessenger.hh (modified) (1 diff)
• PrimaryGeneratorAction.hh (modified) (1 diff)
• PrimaryGeneratorMessenger.hh (modified) (1 diff)
• RunAction.hh (modified) (5 diffs)
• StepMax.hh (modified) (1 diff)
• StepMaxMessenger.hh (modified) (1 diff)
• SteppingAction.hh (modified) (1 diff)
• SteppingVerbose.hh (modified) (1 diff)
• TrackingAction.hh (modified) (1 diff)
• c2_function.hh (modified) (50 diffs)
• c2_function.icc (modified) (32 diffs)

trunk/examples/extended/electromagnetic/TestEm7/include/DetectorConstruction.hh

@@ -25 +25 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/DetectorMessenger.hh

@@ -25 +25 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/EventAction.hh

@@ -25 +25 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/EventActionMessenger.hh

@@ -25 +25 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/G4ScreenedNuclearRecoil.hh

@@ -25 +25 @@
 //
 //
-// $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 +45 @@
 //
 // 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 +55 @@
 #include "globals.hh"
 #include "G4VDiscreteProcess.hh"
+#include "G4ParticleChange.hh"
 #include "c2_function.hh"
 
@@ -59 +61 @@
 #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 +76 @@
 {
 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 +119 @@
         
         // 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 +127 @@
         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 +139 @@
         G4ScreenedCoulombCrossSection *crossSection;
         G4double a1, a2, sinTheta, cosTheta, sinZeta, cosZeta, eRecoil;
-        G4ParticleDefinition *recoilIon;        } G4CoulombKinematicsInfo;
+        G4ParticleDefinition *recoilIon;        
+        const G4Material *targetMaterial;
+} G4CoulombKinematicsInfo;
 
 class G4ScreenedCollisionStage {
@@ -146 +153 @@
 
 public:
-        G4ScreenedCoulombClassicalKinematics() { }
+        G4ScreenedCoulombClassicalKinematics();
         virtual void DoCollisionStep(class G4ScreenedNuclearRecoil *master,
                 const class G4Track& aTrack, const class G4Step& aStep);
@@ -153 +160 @@
                 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 +178 @@
 };
 
+/**
+        \brief A process which handles screened Coulomb collisions between nuclei
+ 
+*/
+
 class G4ScreenedNuclearRecoil : public G4ScreenedCoulombCrossSectionInfo, public G4VDiscreteProcess
 {
@@ -171 +189 @@
         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 +312 @@
         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 +337 @@
         
         std::map<G4int, G4ScreenedCoulombCrossSection*> crossSectionHandlers;           
-        std::map<G4int, c2_function<G4double>*> meanFreePathTables;
         
         G4bool validCollision;
         G4CoulombKinematicsInfo kinematics;
-
+        const G4VNIELPartition *NIELPartitionFunction;
 };
 
@@ -268 +361 @@
         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 +372 @@
 };
 
-
+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

trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmLivermore.hh

@@ -26 +26 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmPenelope.hh

@@ -26 +26 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmStandard.hh

@@ -26 +26 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmStandardNR.hh

@@ -25 +25 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/PhysListEmStandardSS.hh

@@ -25 +25 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/PhysicsList.hh

@@ -25 +25 @@
 //
 //
-// $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 +39 @@
 
 #include "G4VModularPhysicsList.hh"
+#include "G4EmConfigurator.hh"
 #include "globals.hh"
 
@@ -50 +51 @@
 {
 public:
+
   PhysicsList();
   virtual ~PhysicsList();
@@ -67 +69 @@
 
 private:
+
+  G4EmConfigurator em_config; 
+
   G4double cutForGamma;
   G4double cutForElectron;
@@ -77 +82 @@
   G4String                             emName;
   G4VPhysicsConstructor*               emPhysicsList;
+  G4VPhysicsConstructor*               decPhysicsList;
   std::vector<G4VPhysicsConstructor*>  hadronPhys;
     

trunk/examples/extended/electromagnetic/TestEm7/include/PhysicsListMessenger.hh

@@ -25 +25 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/PrimaryGeneratorAction.hh

@@ -25 +25 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/PrimaryGeneratorMessenger.hh

@@ -25 +25 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/RunAction.hh

@@ -24 +24 @@
 // ********************************************************************
 //
-// $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 +40 @@
 class PrimaryGeneratorAction;
 class G4Run;
-
-namespace AIDA {
- class IAnalysisFactory;
- class ITree;
- class IHistogram1D;
-} 
+class Histo;
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -52 +47 @@
 {
 public:
-  RunAction(DetectorConstruction*, PhysicsList*,PrimaryGeneratorAction*);
+  RunAction(DetectorConstruction*, PhysicsList*,
+            PrimaryGeneratorAction*);
   virtual ~RunAction();
 
@@ -58 +54 @@
   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 +79 @@
   G4double                edeptot, eniel;
   G4int                   nPrimarySteps;
-           
-  AIDA::IAnalysisFactory* af;  
-  AIDA::ITree*            tree;
-  AIDA::IHistogram1D*     histo[1];        
+  G4int                   nRange;
+
+  Histo*                  histo;
 };
 

trunk/examples/extended/electromagnetic/TestEm7/include/StepMax.hh

@@ -25 +25 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/StepMaxMessenger.hh

@@ -25 +25 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/SteppingAction.hh

@@ -25 +25 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/SteppingVerbose.hh

@@ -25 +25 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/TrackingAction.hh

@@ -25 +25 @@
 //
 // $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......

trunk/examples/extended/electromagnetic/TestEm7/include/c2_function.hh

@@ +1 @@
+//
+// ********************************************************************
+// * 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 +33 @@
  *  \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 +74 @@
 
 // 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 +128 @@
     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 +138 @@
     /// \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 +147 @@
 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 +173 @@
         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 +200 @@
     /// \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 +222 @@
     /// 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 +236 @@
         /// \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 +245 @@
     /// 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 +300 @@
         /// \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 +313 @@
         /// \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 +325 @@
         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 +334 @@
         /// \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 +351 @@
         /// \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 +357 @@
         /// \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 +364 @@
         /// \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 +441 @@
                 }
         
+        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 +505 @@
     /// 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 +771 @@
 /// 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 +798 @@
         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 +902 @@
         /// \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 +921 @@
 
 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 +958 @@
 
 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 +966 @@
 
 /// \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 +977 @@
         /// \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 +1002 @@
 
 /// \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 +1129 @@
 };
 
-/// \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 +1144 @@
 };
 
+/// \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 +1359 @@
     (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 +1370 @@
     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 +1396 @@
     /// \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 +1510 @@
         // 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 +1524 @@
     /// \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 +1658 @@
     /// \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 +1692 @@
                 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 +1753 @@
         /// \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 +1787 @@
     /// \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 +1800 @@
                 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 +1814 @@
  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 +1823 @@
     /// \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 +1837 @@
                 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 +1869 @@
                 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 +1885 @@
  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 +1906 @@
     ///  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 +1951 @@
 */
 
-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 +1965 @@
 
 /**
-  \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 +2018 @@
 /// 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 +2109 @@
         /// 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;
 };

trunk/examples/extended/electromagnetic/TestEm7/include/c2_function.icc

@@ +1 @@
+//
+// ********************************************************************
+// * 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 +33 @@
  *  \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 +48 @@
 
 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 +59 @@
         // 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 +72 @@
         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 +143 @@
            }
            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 +174 @@
 
 // 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 +204 @@
                 
                 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 +303 @@
                 }
                 
-                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 +364 @@
 }
 
-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 +410 @@
                 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 +447 @@
 }
 
-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 +454 @@
         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 +523 @@
         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 +534 @@
         
         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 +544 @@
                 
                 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 +616 @@
 
         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 +636 @@
         }
         
-        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 +650 @@
                 "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 +812 @@
         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 +918 @@
         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 +943 @@
                 }
         } 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 +968 @@
         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 +1004 @@
 }
 
-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 +1023 @@
 }
 
-// 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 +1028 @@
 // 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 +1056 @@
         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 +1099 @@
 }
 
-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 +1162 @@
                 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 +1173 @@
 }
 
-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 +1207 @@
 }
 
-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 +1219 @@
                         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);
+}