---

func.h

Go to the documentation of this file.

#ifndef __func_h
#define __func_h

#include <math.h>
#include <string>
#include <iostream>
#include <sstream>
#include <mak/infinity.h>
#include "punnets_base.h"

#define USE_DYNAMIC

namespace punnets_common {


class message_base
{
public:
    message_base() { }
    virtual ~message_base() { }

    virtual const char *getMessageId() const = 0;
};  


class func_base
{
protected:
    virtual void valueChange() { }
    
public:
    func_base() { }
    virtual ~func_base() { }

    virtual bool shouldDelete(ntime_t /*current*/) { return false; }
    virtual bool processMessage(ntime_t, const message_base &) { return false; }

    virtual void setLambda(real /*new_lambda*/) { }
    virtual void setZeroPoint(real /*new_zeropoint*/) { }

    virtual real getMaxGradient( ntime_t t ) const = 0;
    virtual ntime_t getNextIncontinuity( ntime_t /*from*/ ) const { return mak::Infinity; };

    virtual real getValue( ntime_t t ) const = 0;
    virtual real get1stDeriv( ntime_t t ) const = 0;
    virtual real get2ndDeriv( ntime_t t ) const = 0;

    virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const = 0;
    virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const = 0;
    virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const = 0;

    virtual func_base * clone() = 0;
    virtual std::string getDescription() = 0;
};  


class func_deriveq_base : public func_base
{
protected:
    real lambda;
    ntime_t zerop;

    virtual void zeropChange() { }
    
public:
    func_deriveq_base() : func_base(), lambda(1.0), zerop(-10000000.0) { }
    virtual ~func_deriveq_base() { }

    virtual void setLambda(real new_lambda) { assert( new_lambda > 0.0 ); lambda = new_lambda; valueChange(); }
    virtual void setZeroPoint(real new_zeropoint) { zerop = new_zeropoint; zeropChange(); }

    class message_set_lambda : public message_base
    {
        real lambda;
public: 
        message_set_lambda(real newlambda) : message_base() { lambda = newlambda; }
        static const char * const messageId;
        virtual const char *getMessageId() const { return messageId; }
        real getLambda() const { return lambda; }
    };
    class message_set_zero_point : public message_base
    {
        real zerop;
public: 
        message_set_zero_point(real newzerop) : message_base() { zerop = newzerop; }
        static const char * const messageId;
        virtual const char *getMessageId() const { return messageId; }
        real getZeroPoint() const { return zerop; }
    };

    virtual bool processMessage(ntime_t t, const message_base &m) { 
#ifdef USE_DYNAMIC
        if( const message_set_lambda *mm = dynamic_cast<const message_set_lambda *>(&m) ) {
            setLambda( mm->getLambda() );
            return true;
        }
#else
        if( m.getMessageId() == message_set_lambda::messageId ) {
            setLambda( dynamic_cast<const message_set_lambda &>(m).getLambda() );
            return true;
        }
#endif
#ifdef USE_DYNAMIC
        if( const message_set_zero_point *mm = dynamic_cast<const message_set_zero_point *>(&m) ) {
            setZeroPoint( mm->getZeroPoint() );
            return true;
        }
#else
        if( m.getMessageId() == message_set_zero_point::messageId ) {
            setZeroPoint( dynamic_cast<const message_set_zero_point &>(m).getZeroPoint() );
            return true;
        }
#endif
        return func_base::processMessage(t,m); 
    }
};  



class func_const : public func_base
{
    real c;
public:
    func_const(real ic) : func_base(), c(ic) { }            

    virtual real getMaxGradient( ntime_t /*t*/ ) const { return 0; };

    virtual real getValue( ntime_t /*t*/ ) const 
                { return c; }

    virtual real get1stDeriv( ntime_t /*t*/ ) const
                { return 0; }
    virtual real get2ndDeriv( ntime_t /*t*/ ) const
                { return 0; }

    virtual void getValueDomain( ntime_t /*t*/, real &upslope, real &ceil, real &downslope, real &floor ) const 
                { upslope = downslope = 0; ceil = floor = c;  }

    virtual void get1stDerivDomain( ntime_t /*t*/, real &upslope, real &ceil, real &downslope, real &floor ) const 
                { upslope = downslope = 0; ceil = floor = 0;  }

    virtual void get2ndDerivDomain( ntime_t /*t*/, real &upslope, real &ceil, real &downslope, real &floor ) const 
                { upslope = downslope = 0; ceil = floor = 0;  }

    virtual func_const * clone() { return new func_const(*this); }
    virtual std::string getDescription() { std::ostringstream oss; oss << "func_const c=" << c; return oss.str(); }
};



class func_const_int : public func_deriveq_base
{
    real c;

    real limit;

    virtual void valueChange() { limit = c / lambda; }
public:
    func_const_int(real ic) : func_deriveq_base(), c(ic) { valueChange(); }         

    virtual real getMaxGradient( ntime_t /*t*/ ) const { return limit < 0 ? 0 : c; };

    virtual real getValue( ntime_t t ) const 
                { return limit * (1 - exp(- lambda * (t - zerop))); }

    virtual real get1stDeriv( ntime_t t ) const
                { return limit * lambda * exp(- lambda * (t - zerop)); }
    virtual real get2ndDeriv( ntime_t t ) const
                { return - limit * lambda * lambda * exp(- lambda * (t - zerop)); }
    virtual real get3rdDeriv( ntime_t t ) const
                { return limit * lambda * lambda * lambda * exp(- lambda * (t - zerop)); }

    virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
    virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
    virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;

    virtual func_const_int * clone() { return new func_const_int(*this); }
    virtual std::string getDescription() { std::ostringstream oss; oss << "func_const_int c=" << c << ", lambda=" << lambda << ", zerop=" << zerop; return oss.str(); }
};



class func_step : public func_base
{
    real r;
    ntime_t t0;
public:
    func_step(real ir, ntime_t it0) : func_base(), r(ir), t0(it0) { }           

    virtual real getMaxGradient( ntime_t /*t*/ ) const { return 0; };

    virtual real getValue( ntime_t t ) const 
                { return (t0 < t ? r : 0.0); }

    virtual real get1stDeriv( ntime_t /*t*/ ) const
                { return 0; }
    virtual real get2ndDeriv( ntime_t /*t*/ ) const
                { return 0; }

    virtual ntime_t getNextIncontinuity( ntime_t from ) const { return from < t0 ? t0 : mak::Infinity; };

    virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
                { upslope = downslope = 0; ceil = floor = (t >= t0 ? r : 0);  }

    virtual void get1stDerivDomain( ntime_t /*t*/, real &upslope, real &ceil, real &downslope, real &floor ) const 
                { upslope = downslope = 0; ceil = floor = 0;  }

    virtual void get2ndDerivDomain( ntime_t /*t*/, real &upslope, real &ceil, real &downslope, real &floor ) const 
                { upslope = downslope = 0; ceil = floor = 0;  }

    virtual func_step * clone() { return new func_step(*this); }
    virtual std::string getDescription() { std::ostringstream oss; oss << "func_step r=" << r << ", t0=" << t0; return oss.str(); }
};


class func_delta_int : public func_deriveq_base
{
    real r;
    ntime_t t0;

    virtual void valueChange() { }

public:
    func_delta_int(real ir, ntime_t it0) : func_deriveq_base(), r(ir), t0(it0) { valueChange(); }           

    virtual real getMaxGradient( ntime_t /*t*/ ) const { return r<0 ? - lambda * r : 0 ; };

    void setParam(real ir, ntime_t it0) { r = ir; t0 = it0; }

    virtual real getValue( ntime_t t ) const 
                { 
//                  std::cout << "zerop=" << zerop << ", t0=" << t0 << ",t=" << t << std::endl;
                    return (zerop < t0 && t0 <= t) ? r * exp( - lambda * (t - t0) ) : 0.0; }

    virtual real get1stDeriv( ntime_t t ) const
                { return (zerop < t0 && t0 <= t) ? - lambda * r * exp( - lambda * (t - t0) ) : 0.0; }
    virtual real get2ndDeriv( ntime_t t ) const
                { return (zerop < t0 && t0 <= t) ? lambda * lambda * r * exp( - lambda * (t - t0) ) : 0.0; }
    virtual real get3rdDeriv( ntime_t t ) const
                { return (zerop < t0 && t0 <= t) ? - lambda * lambda * lambda * r * exp( - lambda * (t - t0) ) : 0.0; }

    virtual ntime_t getNextIncontinuity( ntime_t from ) const { return from < t0 ? t0 : mak::Infinity; };

    virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
    virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
    virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;

    virtual std::string getDescription() { std::ostringstream oss; oss << "func_delta_int r=" << r << ", t0=" << t0 << ", lambda=" << lambda << ", zerop=" << zerop; return oss.str(); }

    class message_add_pulse : public message_base
    {
        real pulselev;
public: 
        message_add_pulse(real ipl) : message_base() { pulselev = ipl; }
        static const char * const messageId;
        virtual const char *getMessageId() const { return messageId; }
        real getPulseLevel() const { return pulselev; }
    };

    virtual bool processMessage(ntime_t t, const message_base &m) { 
#ifdef USE_DYNAMIC
        if( const message_add_pulse *mm = dynamic_cast<const message_add_pulse *>(&m) ) {
            setParam( getValue(t) + mm->getPulseLevel(), t );
            return true;
        }
#else
        if( m.getMessageId() == message_add_pulse::messageId ) {
            setParam( getValue(t) + dynamic_cast<const message_add_pulse &>(m).getPulseLevel(), t );
            return true;
        }
#endif
        return func_deriveq_base::processMessage(t, m); 
    }

    virtual func_delta_int * clone() { return new func_delta_int(*this); }
};


class func_response : public func_deriveq_base
{
    real r, rinit;

    virtual void valueChange() { }

public:
    func_response(real ir) : func_deriveq_base(), r(ir), rinit(ir) { valueChange(); }           

    virtual real getMaxGradient( ntime_t /*t*/ ) const { return r<0 ? - lambda * r : 0 ; };

    virtual real getValue( ntime_t t ) const 
                { 
                    return r * exp( - lambda * (t - zerop) ) ; }

    virtual real get1stDeriv( ntime_t t ) const
                { return - lambda * r * exp( - lambda * (t - zerop) ); }
    virtual real get2ndDeriv( ntime_t t ) const
                { return lambda * lambda * r * exp( - lambda * (t - zerop) ); }
    virtual real get3rdDeriv( ntime_t t ) const
                { return - lambda * lambda * lambda * r * exp( - lambda * (t - zerop) ); }

    virtual ntime_t getNextIncontinuity( ntime_t /*from*/ ) const { return mak::Infinity; };

    virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
    virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
    virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;

    virtual std::string getDescription() { std::ostringstream oss; oss << "func_response r=" << r << ", lambda=" << lambda << ", zerop=" << zerop; return oss.str(); }

    virtual func_response * clone() { return new func_response(*this); }

    virtual void setZeroPoint(real new_zeropoint) { 
            r = rinit + r * exp( - lambda * (new_zeropoint - zerop) );
                zerop = new_zeropoint; zeropChange(); }
};




class func_sine : public func_base
{
    real r;
    real omega;
    real theta;

    static const real alpha = 1.311;
    static const real beta  = 0.375867;

    void valueChange() { }

public:
    func_sine(real ir, real iomega, real itheta) : func_base(), r(ir), omega(iomega), theta(itheta) { valueChange(); }          

    virtual real getMaxGradient( ntime_t /*t*/ ) const { return omega * r; };

    virtual real getValue( ntime_t t ) const
                { return r * sin(omega * t + theta); }

    virtual real get1stDeriv( ntime_t t ) const
                { return omega * r * cos(omega * t + theta); }

    virtual real get2ndDeriv( ntime_t t ) const
                { return - omega * omega * r * sin(omega * t + theta); }

    virtual real get3rdDeriv( ntime_t t ) const
                { return - omega * omega * omega * r * cos(omega * t + theta); }

    virtual ntime_t getNextIncontinuity( ntime_t /*from*/ ) const { return mak::Infinity; };

    virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
    virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
    virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;

    virtual func_sine * clone() { return new func_sine(*this); }
    virtual std::string getDescription() { std::ostringstream oss; oss << "func_sine r=" << r << ", omega=" << omega << ", theta=" << theta; return oss.str(); }
};



class func_sineshot : public func_base
{
    real r;
    real omega;
    real t0;

    real r_div_2;       
    real omega_2;       
    real duration;      

    static const real alpha = 1.311;
    static const real beta  = 0.375867;

    void valueChange() { duration = M_PI / omega; omega_2 = 2 * omega; r_div_2 = r * 0.5; }

public:
    func_sineshot(real ir, real iomega, real it0) : func_base(), r(ir), omega(iomega), t0(it0) { valueChange(); }           

    virtual real getMaxGradient( ntime_t /*t*/ ) const { return fabs(r_div_2) * omega_2; };

    virtual bool shouldDelete(ntime_t current) { return (current >= t0 + duration); }

    virtual real getValue( ntime_t t ) const
                { return (t0 <= t && t <= t0 + duration ? r_div_2 * (1-cos(omega_2 * (t-t0) )) : 0.0); }

    virtual real get1stDeriv( ntime_t t ) const
                { return (t0 <= t && t <= t0 + duration ? r_div_2 * omega_2 * sin(omega_2 * (t-t0) ) : 0.0); }

    virtual real get2ndDeriv( ntime_t t ) const
                { return (t0 <= t && t <= t0 + duration ? r_div_2 * omega_2 * omega_2 * cos(omega_2 * (t-t0) ) : 0.0); }

    virtual ntime_t getNextIncontinuity( ntime_t from ) const { return from < t0 ? t0 : mak::Infinity; };

    virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
    virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
    virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
    virtual func_sineshot * clone() { return new func_sineshot(*this); }
    virtual std::string getDescription() { std::ostringstream oss; oss << "func_sineshot r=" << r << ", omega=" << omega << ", t0=" << t0; return oss.str(); }

    class message_set_t0 : public message_base
    {
public: 
        message_set_t0() : message_base() { }
        static const char * const messageId;
        virtual const char *getMessageId() const { return messageId; }
    };

    virtual bool processMessage(ntime_t t, const message_base &m) { 
#ifdef USE_DYNAMIC
        if( dynamic_cast<const message_set_t0 *>(&m) != NULL ) 
#else
        if( m.getMessageId() == message_set_t0::messageId ) 
#endif
        {
            t0 = t;
            return true;
        }
        return func_base::processMessage(t, m); 
    }
};



class func_sine_int : public func_deriveq_base
{
    real r;
    real omega;
    real theta;

    real r_inv_lsq_osq; 
    real zerop_value;

    void valueChange() { r_inv_lsq_osq = r / (lambda * lambda + omega * omega); }
    void zeropChange() { zerop_value = lambda * sin( theta + omega * zerop ) - omega * cos( theta + omega * zerop ); }

public:
    func_sine_int(real ir, real iomega, real itheta) : func_deriveq_base(), r(ir), omega(iomega), theta(itheta) { valueChange(); }          

    virtual real getMaxGradient( ntime_t /*t*/ ) const { return r_inv_lsq_osq * (omega * (lambda + 1) - lambda * lambda * zerop_value); };

    virtual real getValue(ntime_t t) const
                { return r_inv_lsq_osq * (
                            lambda * sin( theta + omega * t ) - omega * cos( theta + omega * t )
                            + exp( lambda * (zerop - t) ) * zerop_value );
                }
    virtual real get1stDeriv(ntime_t t) const
                { return r_inv_lsq_osq * (
                            omega * (lambda * cos( theta + omega * t ) + omega * sin( theta + omega * t ))
                                - lambda * exp( lambda * (zerop - t) ) * zerop_value );
                }
    virtual real get2ndDeriv(ntime_t t) const
                { return r_inv_lsq_osq * (
                            omega * omega * (- lambda * sin( theta + omega * t ) + omega * cos( theta + omega * t ))
                                + lambda * lambda * exp( lambda * (zerop - t) ) * zerop_value );
                }
    virtual real get3rdDeriv(ntime_t t) const
                { return r_inv_lsq_osq * (
                            omega * omega * omega * (- lambda * cos( theta + omega * t ) - omega * sin( theta + omega * t ))
                                - lambda * lambda * lambda * exp( lambda * (zerop - t) ) * zerop_value );
                }

    virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
                { upslope = fabs(omega); downslope = -upslope; ceil = r; floor = -r;
                  if( getValue(t) > 0 ) {   upslope = get1stDeriv(t); if(   upslope < 0 ) upslope *= -1; }
                                   else { downslope = get1stDeriv(t); if( downslope > 0 ) downslope *= -1; } }

    virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
                { upslope = omega * omega; downslope = -upslope; ceil = fabs(omega) * r; floor = -ceil;
                  if( get1stDeriv(t) > 0 ) {   upslope = get2ndDeriv(t); if(   upslope < 0 ) upslope *= -1; }
                                      else { downslope = get2ndDeriv(t); if( downslope > 0 ) downslope *= -1; } }

    virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
                { upslope = fabs( omega * omega * omega); downslope = -upslope; ceil = omega * omega * r; floor = -ceil;
                  if( get2ndDeriv(t) > 0 ) {   upslope = get3rdDeriv(t); if(   upslope < 0 ) upslope *= -1; }
                                      else { downslope = get3rdDeriv(t); if( downslope > 0 ) downslope *= -1; } }

    virtual func_sine_int * clone() { return new func_sine_int(*this); }
    virtual std::string getDescription() { std::ostringstream oss; oss << "func_sine_int r=" << r << ", omega=" << omega << ", theta=" << theta << ", lambda=" << lambda << ", zerop=" << zerop; return oss.str(); }
};


class func_sineshot_int : public func_deriveq_base
{
    real r;
    real omega;

    real r_div_2;       
    real omega_2;       
    real duration;      

    static const real alpha = 1.311;
    static const real beta  = 0.375867;

    void valueChange() { duration = M_PI / omega; omega_2 = 2 * omega; r_div_2 = r * 0.5; }
    void zeropChange() {  }

public:
    func_sineshot_int(real ir, real iomega) : func_deriveq_base(), r(ir), omega(iomega) { valueChange(); }          

    virtual real getMaxGradient( ntime_t /*t*/ ) const { return fabs(r_div_2) * omega_2; };

    virtual real getValue( ntime_t t ) const
                { return (zerop <= t && t <= zerop + duration ? r_div_2 * (1-cos(omega_2 * (t-zerop) )) : 0.0); }

    virtual real get1stDeriv( ntime_t t ) const
                { return (zerop <= t && t <= zerop + duration ? r_div_2 * omega_2 * sin(omega_2 * (t-zerop) ) : 0.0); }

    virtual real get2ndDeriv( ntime_t t ) const
                { return (zerop <= t && t <= zerop + duration ? r_div_2 * omega_2 * omega_2 * cos(omega_2 * (t-zerop) ) : 0.0); }

    virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
    virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;
    virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const;

    virtual func_sineshot_int * clone() { return new func_sineshot_int(*this); }
    virtual std::string getDescription() { std::ostringstream oss; oss << "func_sineshot_int r=" << r << ", omega=" << omega; return oss.str(); }
};


class func_exp : public func_base
{
    real r;
    real psi;
    real t0;

public:
    func_exp(real ir, real ipsi, real it0) : func_base(), r(ir), psi(ipsi), t0(it0) { if( psi < 0.0 ) throw "func_exp"; }   

    virtual real getMaxGradient( ntime_t t ) const { return r>0 ? 0 : get1stDeriv(t); }

    virtual real getValue( ntime_t t ) const
                { return r * exp( -psi * (t - t0) ); }

    virtual real get1stDeriv( ntime_t t ) const
                { return -psi * r * exp( -psi * (t - t0) ); }

    virtual real get2ndDeriv( ntime_t t ) const
                { return psi * psi * r * exp( -psi * (t - t0) ); }

    virtual real get3rdDeriv( ntime_t t ) const
                { return - psi * psi * psi * r * exp( -psi * (t - t0) ); }

    virtual ntime_t getNextIncontinuity( ntime_t /*from*/ ) const { return mak::Infinity; };

    virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
                { if( r >= 0 ) {   upslope = 0; ceil =  mak::Infinity; downslope = func_exp::get1stDeriv(t); floor = 0.0; } 
                  else         { downslope = 0; floor= -mak::Infinity;   upslope = func_exp::get1stDeriv(t); ceil  = 0.0; } }

    virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
                { if( r <= 0 ) {   upslope = 0; ceil =  mak::Infinity; downslope = func_exp::get2ndDeriv(t); floor = 0.0; } 
                  else         { downslope = 0; floor= -mak::Infinity;   upslope = func_exp::get2ndDeriv(t); ceil  = 0.0; } }

    virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
                { if( r >= 0 ) {   upslope = 0; ceil =  mak::Infinity; downslope = func_exp::get3rdDeriv(t); floor = 0.0; } 
                  else         { downslope = 0; floor= -mak::Infinity;   upslope = func_exp::get3rdDeriv(t); ceil  = 0.0; } }

    virtual func_exp * clone() { return new func_exp(*this); }
    virtual std::string getDescription() { std::ostringstream oss; oss << "func_exp r=" << r << ", psi=" << psi << ", t0=" << t0; return oss.str(); }
};




class func_exp_int : public func_deriveq_base
{
    real r;
    real psi;
    real t0;
    real r_per_lambda_minus_psi;

    void valueChange() { if( psi < 0.0 || lambda < 0.0 ) abort(); r_per_lambda_minus_psi = r / ( lambda - psi ); }

public:
    func_exp_int(real ir, real ipsi, real it0) : func_deriveq_base(), r(ir), psi(ipsi), t0(it0) { valueChange(); }

    virtual real getMaxGradient( ntime_t /*t*/ ) const { throw "not supported"; }

    virtual real getValue(ntime_t t) const
                { return r_per_lambda_minus_psi * ( 
                                exp( - psi * (t - t0) )
                                - exp( lambda * (zerop - t) - psi * (zerop - t0) ) );
                }

    virtual real get1stDeriv(ntime_t t) const
                { return r_per_lambda_minus_psi * ( 
                                - psi * exp( - psi * (t - t0) )
                                + lambda * exp( lambda * (zerop - t) - psi * (zerop - t0) ) );
                }

    virtual real get2ndDeriv(ntime_t t) const
                { return r_per_lambda_minus_psi * ( 
                                psi*psi * exp( - psi * (t - t0) )
                                - lambda*lambda * exp( lambda * (zerop - t) - psi * (zerop - t0) ) );
                }

    virtual real get3rdDeriv(ntime_t t) const
                { return r_per_lambda_minus_psi * ( 
                                -psi*psi*psi * exp( - psi * (t - t0) )
                                +lambda*lambda*lambda * exp( lambda * (zerop - t) - psi * (zerop - t0) ) );
                }

    virtual ntime_t getNextIncontinuity( ntime_t /*from*/ ) const { return mak::Infinity; };

    virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
                { 
                    real value = getValue(t);
                    real deriv = get1stDeriv(t);
                    if( r >= 0 )
                    {
                        if( deriv > 0 )
                        {
                            real log_lambda_per_psi = log(lambda / psi) / (lambda - psi);
                            upslope = deriv; 
                            ceil = r_per_lambda_minus_psi * exp( - psi * (zerop - t0) ) * ( exp( -psi * log_lambda_per_psi ) - exp( - lambda * log_lambda_per_psi) );
                            downslope = - value / log_lambda_per_psi; /* can be more efficient */
                            floor = 0;
                        }
                        else {
                            upslope = 0; ceil = mak::Infinity;
                            downslope = -1; throw "Not supported";
                            floor = 0;
                        }
                    }
                    else /* r < 0 */ {
                        if( deriv < 0 )
                        {
                            real log_lambda_per_psi = log(lambda / psi) / (lambda - psi);
                            downslope = deriv; 
                            floor = r_per_lambda_minus_psi * exp( - psi * (zerop - t0) ) * ( exp( -psi * log_lambda_per_psi ) - exp( - lambda * log_lambda_per_psi) );
                            upslope = - value / log_lambda_per_psi; /* can be more efficient */
                            ceil = 0;
                        }
                        else {
                            downslope = 0; floor = -mak::Infinity;
                            upslope = 1; throw "Not supported";
                            ceil = 0;
                        }
                    }
                }

//  virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
    virtual void get1stDerivDomain( ntime_t , real &, real &, real &, real & ) const 
                { throw "not supported"; }

//  virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
    virtual void get2ndDerivDomain( ntime_t , real &, real &, real &, real & ) const 
                { throw "not supported"; }

    virtual func_exp_int * clone() { return new func_exp_int(*this); }
    virtual std::string getDescription() { std::ostringstream oss; oss << "func_exp_int r=" << r << ",t0=" << t0; return oss.str(); }
};


class func_exp_diff : public func_base
{
    real psi1, psi2;
    real phi, scale, invscale;
    real t0orig;
    real r0, t0;
    real r1, t1;
    real r2, t2;
    real it0, ir1, ir2;
    real r1orig, r2orig;
    real maxval;
    real xi;

    void valueChange();

public:
    func_exp_diff(real ir1, real ipsi1, real ir2, real ipsi2, real it0);

    virtual real getMaxGradient( ntime_t t ) const { return r0<0 ? r0 * xi : (t >= t1 ? 0 : get1stDeriv(t)); }

    virtual real getValue( ntime_t t ) const
                { return r0 * (exp( -psi1 * (t - t0) ) - exp( -psi2 * (t - t0) ) ); }

    virtual real get1stDeriv( ntime_t t ) const
                { return r1 * (exp( -psi1 * (t - t1) ) - exp( -psi2 * (t - t1) ) ); }

    virtual real get2ndDeriv( ntime_t t ) const
                { return r2 * (exp( -psi1 * (t - t2) ) - exp( -psi2 * (t - t2) ) ); }

    virtual real get3rdDeriv( ntime_t t ) const
                { return r2 * (-psi1 * exp( -psi1 * (t - t2) ) + psi2 * exp( -psi2 * (t - t2) ) ); }

    virtual ntime_t getNextIncontinuity( ntime_t /*from*/ ) const { return mak::Infinity; }

    virtual void getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const ;
    virtual void get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const ;
    virtual void get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const ;

    virtual func_exp_diff * clone() { return new func_exp_diff(*this); }
    virtual std::string getDescription() { std::ostringstream oss; oss << "func_exp_diff r0=" << r0 << ", psi1=" << psi1 << ", psi2=" << psi2 << ", t0=" << t0; return oss.str(); }


    class message_add_event_time : public message_base
    {
    public:
        message_add_event_time() { }
        static const char * const messageId;
        virtual const char *getMessageId() const { return messageId; }
    };  

    virtual bool processMessage(ntime_t t, const message_base &m);
};


} // namespace punnets

#endif // __func_h

---