---

func.cpp

Go to the documentation of this file.

#include "func.h"
#include <iostream>

namespace punnets_common {

using namespace mak;


const char * const func_deriveq_base::message_set_lambda::messageId = "message_set_lambda";
const char * const func_deriveq_base::message_set_zero_point::messageId = "message_set_zerop";
const char * const func_delta_int::message_add_pulse::messageId = "message_add_pulse";
const char * const func_sineshot::message_set_t0::messageId = "message_set_t0";
const char * const func_exp_diff::message_add_event_time::messageId = "message_add_event_time";

/***************************************************************************/
/**************** func_const_int *******************************************/
/***************************************************************************/

void func_const_int::getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
    if( c > 0 ) { upslope = func_const_int::get1stDeriv(t); downslope = 0; ceil = limit; floor = 0;  }
    else        { downslope = func_const_int::get1stDeriv(t); upslope = 0; floor = limit; ceil = 0;  }
}

void func_const_int::get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
    if( c > 0 ) { upslope = 0; downslope = func_const_int::get2ndDeriv(t); ceil = limit * lambda; floor = 0;  }
    else        { downslope = 0; upslope = func_const_int::get2ndDeriv(t); floor = limit * lambda; ceil = 0;  }
}

void func_const_int::get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
    if( c > 0 ) { upslope = func_const_int::get3rdDeriv(t); downslope = 0; ceil = 0; floor = -Infinity;  }
    else        { downslope = func_const_int::get3rdDeriv(t); upslope = 0; floor = 0; ceil = Infinity;  }
}

/***************************************************************************/
/**************** func_delta_int *******************************************/
/***************************************************************************/

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

void func_delta_int::get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
    if( zerop >= t0 || t0 > t ) { upslope = downslope = 0; ceil = floor = 0; }  
    else if( r <= 0 ) {   upslope = 0; ceil =  Infinity; downslope = get2ndDeriv(t); floor = 0.0; } 
    else              { downslope = 0; floor= -Infinity;   upslope = get2ndDeriv(t); ceil  = 0.0; } 
}

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

/***************************************************************************/
/**************** func_response *******************************************/
/***************************************************************************/

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

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

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

/***************************************************************************/
/**************** func_sine ************************************************/
/***************************************************************************/

void func_sine::getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
    real phi = omega * t + theta;
    if( phi < 0.0 ) phi -= ::floor(phi / (2*M_PI)) * 2*M_PI;
    phi = fmod(phi, 2*M_PI);
    phi -= M_PI;
    if( phi + 0.5*M_PI >= 0 )
        upslope = omega * r * cos( alpha * cos( beta * (phi + 0.5*M_PI)) + 0.5*M_PI - alpha);
    else
        upslope = omega * r * cos(omega * t + theta);
    phi += M_PI;
    if( phi >= M_PI ) phi -= 2*M_PI;
    if( phi + 0.5*M_PI >= 0 )
        downslope = - omega * r * cos( alpha * cos( beta * (phi + 0.5*M_PI)) + 0.5*M_PI - alpha);
    else
        downslope = omega * r * cos(omega * t + theta);
    ceil = r; floor = -r; 
}

void func_sine::get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
    real phi = omega * t + theta + 0.5*M_PI;
    if( phi < 0.0 ) phi -= ::floor(phi / (2*M_PI)) * 2*M_PI;
    phi = fmod(phi, 2*M_PI);
    phi -= M_PI;
    if( phi + 0.5*M_PI >= 0 )
        upslope = omega * omega * r * cos( alpha * cos( beta * (phi + 0.5*M_PI)) + 0.5*M_PI - alpha);
    else
        upslope = -omega * omega * r * cos(phi);
    phi += M_PI;
    if( phi >= M_PI ) phi -= 2*M_PI;
    if( phi + 0.5*M_PI >= 0 )
        downslope = - omega * omega * r * cos( alpha * cos( beta * (phi + 0.5*M_PI)) + 0.5*M_PI - alpha);
    else
        downslope = omega * omega * r * cos(phi);
    ceil = r*omega; floor = -r*omega; 
}

void func_sine::get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{
    real phi = omega * t + theta + M_PI;
    if( phi < 0.0 ) phi -= ::floor(phi / (2*M_PI)) * 2*M_PI;
    phi = fmod(phi, 2*M_PI);
    phi -= M_PI;
    if( phi + 0.5*M_PI >= 0 )
        upslope = omega * omega * omega * r * cos( alpha * cos( beta * (phi + 0.5*M_PI)) + 0.5*M_PI - alpha);
    else
        upslope = -omega * omega * omega * r * cos(phi);
    phi += M_PI;
    if( phi >= M_PI ) phi -= 2*M_PI;
    if( phi + 0.5*M_PI >= 0 )
        downslope = - omega * omega * omega * r * cos( alpha * cos( beta * (phi + 0.5*M_PI)) + 0.5*M_PI - alpha);
    else
        downslope = omega * omega * omega * r * cos(phi);
    ceil = r*omega*omega; floor = -r*omega*omega; 
}


/***************************************************************************/
/**************** func_sineshot ********************************************/
/***************************************************************************/

void func_sineshot::getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
    upslope = ceil = downslope = floor = 0.0;
    if( t < t0 || t >= t0 + duration )
        return;
    else 
    {
        real phase = omega_2 * (t - t0);
        real phi = phase + M_PI;
        if( phi >= 2*M_PI ) phi -= 2*M_PI;
        ( r > 0 ? ceil : floor ) = r;
        if( phase < 0.5 * M_PI )
        {
            ( r > 0 ? upslope : downslope ) = omega_2 * r_div_2 * cos( alpha * cos( beta * phi ) + 0.5*M_PI - alpha);
            ( r > 0 ? downslope : upslope ) = - omega_2 * r_div_2 * cos( alpha * cos( beta * (phi-M_PI) ) + 0.5*M_PI - alpha);
        }
        else if( phase < M_PI )
        {
            ( r > 0 ? upslope : downslope ) = r_div_2 * omega_2 * sin(phase );
            ( r > 0 ? downslope : upslope ) = - omega_2 * r_div_2 * cos( alpha * cos( beta * (phi-M_PI) ) + 0.5*M_PI - alpha);
        }
        else if( phase < 1.5*M_PI )
        {
            ( r > 0 ? downslope : upslope ) = - omega_2 * r_div_2 * cos( alpha * cos( beta * (phi+M_PI) ) + 0.5*M_PI - alpha);
        }
        else
        {
            ( r > 0 ? downslope : upslope ) = r_div_2 * omega_2 * sin(phase);
        }
    }
}

void func_sineshot::get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
    upslope = ceil = downslope = floor = 0.0;
    if( t < t0 || t >= t0 + duration )
        return;
    else 
    {
        real phase = omega_2 * (t - t0);
        real phi = phase + 1.5*M_PI;
        if( phi >= 2*M_PI ) phi -= 2*M_PI;
        ceil = fabs(r_div_2) * omega_2;
        floor = -ceil;
        if( phase < 0.5 * M_PI )
        {
            ( r > 0 ? upslope : downslope ) = r_div_2 * omega_2 * omega_2 * cos(phase );
            ( r > 0 ? downslope : upslope ) = - omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * (phi - M_PI) ) + 0.5*M_PI - alpha);
        }
        else if( phase < M_PI )
        {
        //  ( r > 0 ? upslope : downslope ) = 0; 
            ( r > 0 ? downslope : upslope ) = - omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * (phi + M_PI) ) + 0.5*M_PI - alpha);
        }
        else if( phase < 1.5*M_PI )
        {
            ( r > 0 ? ceil : floor ) = 0; 
            ( r > 0 ? upslope : downslope ) = (- r_div_2 * omega_2 * sin(phase )) / (duration - (t-t0));
            ( r > 0 ? downslope : upslope ) = r_div_2 * omega_2 * omega_2 * cos(phase );
        }
        else
        {
            ( r > 0 ? ceil : floor ) = 0; 
            ( r > 0 ? upslope : downslope ) = (- r_div_2 * omega_2 * sin(phase )) / (duration - (t-t0));
        }
    }
}

void func_sineshot::get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
    upslope = ceil = downslope = floor = 0.0;
    if( t < t0 || t >= t0 + duration )
        return;
    else 
    {
        real phase = omega_2 * (t - t0);
        real phi = phase;
        ceil = fabs(r_div_2 * omega_2 * omega_2);
        floor = -ceil;
        if( phase < 0.5 * M_PI )
        {
            ( r > 0 ? upslope : downslope ) = omega_2 * omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * phi ) + 0.5*M_PI - alpha);
            ( r > 0 ? downslope : upslope ) = - omega_2 * omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * (phi + M_PI) ) + 0.5*M_PI - alpha);
        }
        else if( phase < M_PI )
        {
            ( r > 0 ? upslope : downslope ) = omega_2 * omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * phi ) + 0.5*M_PI - alpha);
            ( r > 0 ? downslope : upslope ) = - r_div_2 * omega_2 * omega_2 * omega_2 * sin( phase );
        }
        else if( phase < 1.5*M_PI )
        {
            ( r > 0 ? upslope : downslope ) = omega_2 * omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * phi ) + 0.5*M_PI - alpha);
        }
        else
        {
            ( r > 0 ? upslope : downslope ) = - r_div_2 * omega_2 * omega_2 * omega_2 * sin( phase );
            ( r > 0 ? downslope : upslope ) = - (r_div_2 * omega_2 * omega_2 * cos(phase )) / (duration - (t-t0));
            ( r > 0 ? floor : ceil ) = 0; 
        }
    }
}


/***************************************************************************/
/**************** func_sineshot_int ********************************************/
/***************************************************************************/

void func_sineshot_int::getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
    const real t0 = zerop;
    upslope = ceil = downslope = floor = 0.0;
    if( t < t0 || t >= t0 + duration )
        return;
    else 
    {
        real phase = omega_2 * (t - t0);
        real phi = phase + M_PI;
        if( phi >= 2*M_PI ) phi -= 2*M_PI;
        ( r > 0 ? ceil : floor ) = r;
        if( phase < 0.5 * M_PI )
        {
            ( r > 0 ? upslope : downslope ) = omega_2 * r_div_2 * cos( alpha * cos( beta * phi ) + 0.5*M_PI - alpha);
            ( r > 0 ? downslope : upslope ) = - omega_2 * r_div_2 * cos( alpha * cos( beta * (phi-M_PI) ) + 0.5*M_PI - alpha);
        }
        else if( phase < M_PI )
        {
            ( r > 0 ? upslope : downslope ) = r_div_2 * omega_2 * sin(phase );
            ( r > 0 ? downslope : upslope ) = - omega_2 * r_div_2 * cos( alpha * cos( beta * (phi-M_PI) ) + 0.5*M_PI - alpha);
        }
        else if( phase < 1.5*M_PI )
        {
            ( r > 0 ? downslope : upslope ) = - omega_2 * r_div_2 * cos( alpha * cos( beta * (phi+M_PI) ) + 0.5*M_PI - alpha);
        }
        else
        {
            ( r > 0 ? downslope : upslope ) = r_div_2 * omega_2 * sin(phase);
        }
    }
}

void func_sineshot_int::get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
    const real t0 = zerop;
    upslope = ceil = downslope = floor = 0.0;
    if( t < t0 || t >= t0 + duration )
        return;
    else 
    {
        real phase = omega_2 * (t - t0);
        real phi = phase + 1.5*M_PI;
        if( phi >= 2*M_PI ) phi -= 2*M_PI;
        ceil = fabs(r_div_2) * omega_2;
        floor = -ceil;
        if( phase < 0.5 * M_PI )
        {
            ( r > 0 ? upslope : downslope ) = r_div_2 * omega_2 * omega_2 * cos(phase );
            ( r > 0 ? downslope : upslope ) = - omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * (phi - M_PI) ) + 0.5*M_PI - alpha);
        }
        else if( phase < M_PI )
        {
        //  ( r > 0 ? upslope : downslope ) = 0; 
            ( r > 0 ? downslope : upslope ) = - omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * (phi + M_PI) ) + 0.5*M_PI - alpha);
        }
        else if( phase < 1.5*M_PI )
        {
            ( r > 0 ? ceil : floor ) = 0; 
            ( r > 0 ? upslope : downslope ) = (- r_div_2 * omega_2 * sin(phase )) / (duration - (t-t0));
            ( r > 0 ? downslope : upslope ) = r_div_2 * omega_2 * omega_2 * cos(phase );
        }
        else
        {
            ( r > 0 ? ceil : floor ) = 0; 
            ( r > 0 ? upslope : downslope ) = (- r_div_2 * omega_2 * sin(phase )) / (duration - (t-t0));
        }
    }
}

void func_sineshot_int::get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
    const real t0 = zerop;
    upslope = ceil = downslope = floor = 0.0;
    if( t < t0 || t >= t0 + duration )
        return;
    else 
    {
        real phase = omega_2 * (t - t0);
        real phi = phase;
        ceil = fabs(r_div_2 * omega_2 * omega_2);
        floor = -ceil;
        if( phase < 0.5 * M_PI )
        {
            ( r > 0 ? upslope : downslope ) = omega_2 * omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * phi ) + 0.5*M_PI - alpha);
            ( r > 0 ? downslope : upslope ) = - omega_2 * omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * (phi + M_PI) ) + 0.5*M_PI - alpha);
        }
        else if( phase < M_PI )
        {
            ( r > 0 ? upslope : downslope ) = omega_2 * omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * phi ) + 0.5*M_PI - alpha);
            ( r > 0 ? downslope : upslope ) = - r_div_2 * omega_2 * omega_2 * omega_2 * sin( phase );
        }
        else if( phase < 1.5*M_PI )
        {
            ( r > 0 ? upslope : downslope ) = omega_2 * omega_2 * omega_2 * r_div_2 * cos( alpha * cos( beta * phi ) + 0.5*M_PI - alpha);
        }
        else
        {
            ( r > 0 ? upslope : downslope ) = - r_div_2 * omega_2 * omega_2 * omega_2 * sin( phase );
            ( r > 0 ? downslope : upslope ) = - (r_div_2 * omega_2 * omega_2 * cos(phase )) / (duration - (t-t0));
            ( r > 0 ? floor : ceil ) = 0; 
        }
    }
}

/***************************************************************************/
/**************** func_exp_diff ********************************************/
/***************************************************************************/

func_exp_diff::func_exp_diff(real iir1, real ipsi1, real iir2, real ipsi2, real iit0) 
    : func_base(), psi1(ipsi1), psi2(ipsi2), t0orig(iit0), r1orig(iir1), r2orig(iir2)
{
    it0 = t0orig;
    ir1 = r1orig;
    ir2 = r2orig;
    valueChange();
};

void func_exp_diff::valueChange()
{
    if( psi1 < 0.0 || psi2 <= psi1 || ir1 *ir2 <= 0.0 ) 
    {
        std::cout << "Abort!! bad condition!! " << std::endl;
        std::cout << "psi1 = " << psi1 << ", psi2 = "<< psi2 << ", ir1 = " << ir1 << ", ir2 = " << ir2 << std::endl;
        throw "func_exp_diff"; 
    }

    phi = psi2 / psi1;
    real z = 2 * log( phi ) / (phi - 1);
    scale = psi1 / z;
    invscale = z / psi1;
    maxval = exp(-psi1 * 0.5 * invscale) - exp(-psi2 * 0.5 * invscale);
    xi = - psi1 * exp(-psi1 * invscale) + psi2 * exp(-psi2 * invscale);

    t0 = (log(ir1/ir2))/(psi1-psi2);
    r0 = ir1 * exp( -psi1 * t0 ); 

    t1 = (log(ir1*psi1/(ir2*psi2)))/(psi1-psi2);
    r1 = - ir1 * psi1 * exp( -psi1 * t1 );

    t2 = (log(ir1*psi1*psi1/(ir2*psi2*psi2)))/(psi1-psi2);
    r2 = ir1 * psi1 * psi1 * exp( -psi1 * t2 );

    if( fabs( ir1 * exp( -psi1 * t0 ) - ir2 * exp( -psi2 * t0 ) ) > 1e-5 ) 
        std::cout << "komatta 0! " << ir1 * exp( -psi1 * t0 ) << ", " << ir2 * exp( -psi2 * t0 ) << std::endl << getDescription() << std::endl;
    if( fabs( ir1 * psi1 * exp( -psi1 * t1 ) - ir2 * psi2 * exp( -psi2 * t1 ) ) > 1e-5 )
        std::cout << "komatta 1! " << ir1 * psi1 * exp( -psi1 * t1 ) << ", " << ir2 * psi2 * exp( -psi2 * t1 ) << std::endl;
    if( fabs( ir1 * psi1 * psi1 * exp( -psi1 * t2 ) - ir2 * psi2 * psi2 * exp( -psi2 * t2 ) ) > 1e-5 )
        std::cout << "komatta 2! " << ir1 * psi1 * psi1 * exp( -psi1 * t2 ) << ", " << ir2 * psi2 * psi2 * exp( -psi2 * t2 ) << std::endl;
//  cout << "scale = "<< scale << ", maxval = "<< maxval << ", xi= "<< xi << endl;
//  cout << "t0 = "<< t0 << ", r0 = "<< r0 << endl;
//  cout << "t1 = "<< t1 << ", r1 = "<< r1 << endl;
//  cout << "t2 = "<< t2 << ", r2 = "<< r2 << endl;
    t0 += it0; t1 += it0; t2 += it0;
    phi = 0.3 *(log(phi)/log(10)) * (log(phi)/log(10)) + 2.45;
}   

void func_exp_diff::getValueDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{
    real val = getValue(t);
    real der = get1stDeriv(t);

    if( t < t0 + 0.5 * invscale ) 
    {
        if( t < t0 ) 
        {
            if( r0 > 0 ) { ceil  = maxval * r0; floor = val; upslope   = der; downslope = 0; }
            else         { floor = maxval * r0; ceil  = val; downslope = der; upslope = 0;   }
        }
        else // if( t < t0 + 0.5 * invscale ) 
        {
            if( r0 > 0 ) { ceil  = maxval * r0; floor = 0; upslope   = der; downslope = r0 * xi * (1.0 - pow(1.0 - (t-t0)*scale, phi)); }
            else         { floor = maxval * r0; ceil  = 0; downslope = der; upslope   = r0 * xi * (1.0 - pow(1.0 - (t-t0)*scale, phi)); }
        }
    }
    else if( t < t0 + 1.0 * invscale ) 
    {
        if( r0 > 0 ) { ceil  = val; floor = 0; upslope   = 0; downslope = r0 * xi * (1.0 - pow(1.0 - (t-t0)*scale, phi)); }
        else         { floor = val; ceil  = 0; downslope = 0; upslope   = r0 * xi * (1.0 - pow(1.0 - (t-t0)*scale, phi)); }
    }
    else
    {
        if( r0 > 0 ) { ceil  = val; floor = 0; upslope   = 0; downslope = der; }
        else         { floor = val; ceil  = 0; downslope = 0; upslope   = der; }
    }
}

void func_exp_diff::get1stDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{ 
    real val = get1stDeriv(t);
    real der = get2ndDeriv(t);

    if( t < t1 + 0.5 * invscale ) 
    {
        if( t < t1 ) 
        {
            if( r1 > 0 ) { ceil  = maxval * r1; floor = val; upslope   = der; downslope = 0; }
            else         { floor = maxval * r1; ceil  = val; downslope = der; upslope = 0;   }
        }
        else // if( t < t1 + 0.5 * invscale ) 
        {
            if( r1 > 0 ) { ceil  = maxval * r1; floor = 0; upslope   = der; downslope = r1 * xi * (1.0 - pow(1.0 - (t-t1)*scale, phi)); }
            else         { floor = maxval * r1; ceil  = 0; downslope = der; upslope   = r1 * xi * (1.0 - pow(1.0 - (t-t1)*scale, phi)); }
        }
    }
    else if( t < t1 + 1.0 * invscale ) 
    {
        if( r1 > 0 ) { ceil  = val; floor = 0; upslope   = 0; downslope = r1 * xi * (1.0 - pow(1.0 - (t-t1)*scale, phi)); }
        else         { floor = val; ceil  = 0; downslope = 0; upslope   = r1 * xi * (1.0 - pow(1.0 - (t-t1)*scale, phi));   }
    }
    else
    {
        if( r1 > 0 ) { ceil  = val; floor = 0; upslope   = 0; downslope = der; }
        else         { floor = val; ceil  = 0; downslope = 0; upslope   = der; }
    }
}

void func_exp_diff::get2ndDerivDomain( ntime_t t, real &upslope, real &ceil, real &downslope, real &floor ) const 
{
    real val = get2ndDeriv(t);
    real der = get3rdDeriv(t);

    if( t < t2 + 0.5 * invscale ) 
    {
        if( t < t2 ) 
        {
            if( r2 > 0 ) { ceil  = maxval * r2; floor = val; upslope   = der; downslope = 0; }
            else         { floor = maxval * r2; ceil  = val; downslope = der; upslope = 0;   }
        }
        else // if( t < t2 + 0.5 * invscale ) 
        {
            if( r2 > 0 ) { ceil  = maxval * r2; floor = 0; upslope   = der; downslope = r2 * xi * (1.0 - pow(1.0 - (t-t2)*scale, phi)); }
            else         { floor = maxval * r2; ceil  = 0; downslope = der; upslope   = r2 * xi * (1.0 - pow(1.0 - (t-t2)*scale, phi)); }
        }
    }
    else if( t < t2 + 1.0 * invscale ) 
    {
        if( r2 > 0 ) { ceil  = val; floor = 0; upslope   = 0; downslope = r2 * xi * (1.0 - pow(1.0 - (t-t2)*scale, phi)); }
        else         { floor = val; ceil  = 0; downslope = 0; upslope   = r2 * xi * (1.0 - pow(1.0 - (t-t2)*scale, phi));   }
    }
    else
    {
        if( r2 > 0 ) { ceil  = val; floor = 0; upslope   = 0; downslope = der; }
        else         { floor = val; ceil  = 0; downslope = 0; upslope   = der; }
    }
}

bool func_exp_diff::processMessage(ntime_t t, const message_base &m) 
{ 
#ifdef USE_DYNAMIC
    if( dynamic_cast<const message_add_event_time *>(&m) != NULL ) 
#else
    if( m.getMessageId() == message_add_event_time::messageId ) 
#endif
    {
        ir1 = r1orig + ir1 * exp( - psi1 * (t-it0) );
        ir2 = r2orig + ir2 * exp( - psi2 * (t-it0) );
        it0 = t;
        valueChange();
        return true;
    }
    return func_base::processMessage(t, m); 
}


} // namespace punnets

---