diff --git a/SimulationRuntime/c/simulation/math-support/Makefile.in b/SimulationRuntime/c/simulation/math-support/Makefile.in index bab2bddb4df..68b4b88f85c 100644 --- a/SimulationRuntime/c/simulation/math-support/Makefile.in +++ b/SimulationRuntime/c/simulation/math-support/Makefile.in @@ -14,7 +14,7 @@ FFLAGS = -O -fexceptions OBJS = bigden.o enorm.o hybrd1.o nelmead.o qform.o r1updt.o \ biglag.o dogleg.o fdjac1.o hybrj.o newuoa.o qrfac.o trsapp.o \ dpmpar.o hybrd.o lsame.o newuob.o r1mpyq.o update.o \ -simulation_events.o simulation_delay.o simulation_init.o dgesv_aux.o \ +simulation_events.o ringbuffer.o simulation_delay.o simulation_init.o dgesv_aux.o \ HFILES = blaswrap.h \ matrix.h \ diff --git a/SimulationRuntime/c/simulation/math-support/ringbuffer.c b/SimulationRuntime/c/simulation/math-support/ringbuffer.c new file mode 100644 index 00000000000..c76e18a21d1 --- /dev/null +++ b/SimulationRuntime/c/simulation/math-support/ringbuffer.c @@ -0,0 +1,94 @@ +/* + * This file is part of OpenModelica. + * + * Copyright (c) 1998-2008, Linköpings University, + * Department of Computer and Information Science, + * SE-58183 Linköping, Sweden. + * + * All rights reserved. + * + * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THIS OSMC PUBLIC + * LICENSE (OSMC-PL). ANY USE, REPRODUCTION OR DISTRIBUTION OF + * THIS PROGRAM CONSTITUTES RECIPIENT'S ACCEPTANCE OF THE OSMC + * PUBLIC LICENSE. + * + * The OpenModelica software and the Open Source Modelica + * Consortium (OSMC) Public License (OSMC-PL) are obtained + * from Linköpings University, either from the above address, + * from the URL: http://www.ida.liu.se/projects/OpenModelica + * and in the OpenModelica distribution. + * + * This program is distributed WITHOUT ANY WARRANTY; without + * even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH + * IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS + * OF OSMC-PL. + * + * See the full OSMC Public License conditions for more details. + * + */ + +#include "ringbuffer.h" + +#include +#include +#include + +void allocRingBuffer(RINGBUFFER* rb, int sz, int item_size) +{ + rb->first_element = 0; + rb->num_element = 0; + rb->buf_size = sz > 0 ? sz : 1; + rb->item_size = item_size; + rb->buffer = calloc(rb->buf_size, rb->item_size); + assert(rb->buffer); +} + +void freeRingBuffer(RINGBUFFER* rb) +{ + free(rb->buffer); +} + +void* getRingData(RINGBUFFER* rb, int nIndex) +{ + assert(nIndex < rb->num_element); + assert(0 <= nIndex); + return ((char*)rb->buffer)+(((rb->first_element+nIndex)%rb->buf_size)*rb->item_size); +} + +void expandRingBuffer(RINGBUFFER* rb) +{ + int i; + void* temp = calloc(2*rb->buf_size, rb->item_size); + assert(temp); + + for(i=0; inum_element; i++) + memcpy(((char*)temp)+(i*rb->item_size), getRingData(rb, i), rb->item_size); + + free(rb->buffer); + rb->buffer = temp; + rb->buf_size *= 2; + rb->first_element = 0; +} + +void appendRingData(RINGBUFFER* rb, void* value) +{ + if(rb->buf_size < rb->num_element+1) + expandRingBuffer(rb); + + memcpy(((char*)rb->buffer)+(((rb->first_element+rb->num_element)%rb->buf_size)*rb->item_size), value, rb->item_size); + ++rb->num_element; +} + +void dequeueNFirstRingDatas(RINGBUFFER* rb, int n) +{ + assert(n <= rb->num_element); + assert(0 <= n); + rb->first_element = (rb->first_element+n)%rb->buf_size; + rb->num_element -= n; +} + +int ringBufferLength(RINGBUFFER* rb) +{ + return rb->num_element; +} diff --git a/SimulationRuntime/c/simulation/math-support/ringbuffer.h b/SimulationRuntime/c/simulation/math-support/ringbuffer.h index 9c3a3b4a1d4..86a1f005828 100644 --- a/SimulationRuntime/c/simulation/math-support/ringbuffer.h +++ b/SimulationRuntime/c/simulation/math-support/ringbuffer.h @@ -31,70 +31,37 @@ #ifndef RINGBUFFER_H #define RINGBUFFER_H -template - /* * This is an expanding ring buffer. * When it gets full, it doubles in size. * It's basically a queue which has get(ix) instead of get_first()/delete_first(). */ -class ringbuffer { -private: - - T* buffer; - int first_element; - int num_element; - int buf_size; - - T& fast_get(const int nIndex) { - return buffer[(first_element+nIndex)%buf_size]; - } - - void expand_buffer() { - T* nb = new T[buf_size*2]; - for (int i=0; i +#include "simulation_runtime.h" +#include "ringbuffer.h" #include "assert.h" #include "string.h" -#include "simulation_runtime.h" - +#include #include -#include "ringbuffer.h" double tStart = 0; typedef struct _TimeAndValue { - double time; - double value; + double time; + double value; } t_TimeAndValue; typedef struct _ExpressionDelayBuffer { - long currentIndex; - long maxExpressionBuffer; - t_TimeAndValue *expressionDelayBuffer; + long currentIndex; + long maxExpressionBuffer; + t_TimeAndValue *expressionDelayBuffer; } t_ExpressionDelayBuffer; // the delayStructure looks like a matrix (rows = expressionNumber+currentColumnIndex, columns={time, value}) -typedef ringbuffer t_buffer; -t_buffer **delayStructure; +RINGBUFFER **delayStructure; extern "C" { -extern const int numDelayExpressionIndex; + extern const int numDelayExpressionIndex; } + void initDelay(double startTime) { - // get the start time of the simulation: time.start. - tStart = startTime; - // allocate the memory for rows - delayStructure = new t_buffer*[numDelayExpressionIndex]; - for (int i=0; i= LOG_SOLVER) - { - fprintf(stdout, "initDelay called with startTime = %f\n", startTime); fflush(NULL); - } + int i; + + // get the start time of the simulation: time.start. + tStart = startTime; + + // allocate the memory for rows + delayStructure = new RINGBUFFER*[numDelayExpressionIndex]; + + for(i=0; i= LOG_SOLVER) + { + fprintf(stdout, "initDelay called with startTime = %f\n", startTime); + fflush(NULL); + } } void deinitDelay() { - delete[] delayStructure; + int i; + + for(i=0; i time) - end = i; - else - start = i; - } while (t != time && end > start + 1); - return start; + int start = 0; + int end = ringBufferLength(&delayStruct); + double t; + do + { + int i = (start + end) / 2; + t = ((t_TimeAndValue*)getRingData(&delayStruct, i))->time; + if(t > time) + end = i; + else + start = i; + }while (t != time && end > start + 1); + return start; } void storeDelayedExpression(int exprNumber, double exprValue) { - double time = globalData->timeValue; + t_TimeAndValue tpl; + double time = globalData->timeValue; + + if(exprNumber < 0) + { + // fprintf(stderr, "storeDelayedExpression: Invalid expression number %d.\n", exprNumber); + return; + } - if (exprNumber < 0) { - // fprintf(stderr, "storeDelayedExpression: Invalid expression number %d.\n", exprNumber); - return; - } + if(time < tStart) + { + // fprintf(stderr, "storeDelayedExpression: Time is smaller than starting time. Value ignored.\n"); + return; + } - if (time < tStart) { - // fprintf(stderr, "storeDelayedExpression: Time is smaller than starting time. Value ignored.\n"); - return; - } + // fprintf(stderr, "storeDelayed[%d] %g:%g\n", exprNumber, time, exprValue); - // fprintf(stderr, "storeDelayed[%d] %g:%g\n", exprNumber, time, exprValue); - - // Allocate more space for expressions - assert(exprNumber < numDelayExpressionIndex); + // Allocate more space for expressions + assert(exprNumber < numDelayExpressionIndex); - t_buffer* delayStruct = delayStructure[exprNumber]; - t_TimeAndValue tpl; - tpl.time = time; - tpl.value = exprValue; - delayStruct->append(tpl); + tpl.time = time; + tpl.value = exprValue; + appendRingData(delayStructure[exprNumber], &tpl); } double delayImpl(int exprNumber, double exprValue, double time, double delayTime, double delayMax /* Unused */) { - // fprintf(stderr, "delayImpl: exprNumber = %d, exprValue = %lf, time = %lf, delayTime = %lf\n", exprNumber, exprValue, time, delayTime); - - // Check for errors - assert(exprNumber >= 0); - assert(exprNumber < numDelayExpressionIndex); - - if (time <= tStart) { - // fprintf(stderr, "delayImpl: Entered at time < starting time: %g.\n", exprValue); - return exprValue; - } - - if (delayTime < 0.0) { - throw TerminateSimulationException(globalData->timeValue, - std::string("Negative delay requested.\n")); - } - - t_buffer* delayStruct = delayStructure[exprNumber]; - int length = delayStruct->length(); - - if (length == 0) { - // This occurs in the initialization phase - // fprintf(stderr, "delayImpl: Missing initial value, using argument value %g instead.\n", exprValue); - return exprValue; - } - - // Returns: expr(time?delayTime) for time>time.start + delayTime and - // expr(time.start) for time <= time.start + delayTime. - // The arguments, i.e., expr, delayTime and delayMax, need to be subtypes of Real. - // DelayMax needs to be additionally a parameter expression. - // The following relation shall hold: 0 <= delayTime <= delayMax, - // otherwise an error occurs. If delayMax is not supplied in the argument list, - // delayTime need to be a parameter expression. See also Section 3.7.2.1. - // For non-scalar arguments the function is vectorized according to Section 10.6.12. - - if (time <= tStart + delayTime) { - double res = (*delayStruct)[0].value; - // fprintf(stderr, "findTime: time <= tStart + delayTime: [%d] = %g\n",exprNumber,res); - return res; - } - else { - // return expr(time-delayTime) - assert(delayTime >= 0.0); - double timeStamp = time - delayTime; - double time0, time1, value0, value1; - - int i; + RINGBUFFER* delayStruct = delayStructure[exprNumber]; + + // fprintf(stderr, "delayImpl: exprNumber = %d, exprValue = %lf, time = %lf, delayTime = %lf\n", exprNumber, exprValue, time, delayTime); - // find the row for the lower limit - if (timeStamp > (*delayStruct)[length - 1].time) { - // delay between the last accepted time step and the current time - time0 = (*delayStruct)[length - 1].time; - value0 = (*delayStruct)[length - 1].value; - time1 = time; - value1 = exprValue; + // Check for errors + assert(0 <= exprNumber); + assert(exprNumber < numDelayExpressionIndex); + + if(time <= tStart) + { + // fprintf(stderr, "delayImpl: Entered at time < starting time: %g.\n", exprValue); + return exprValue; } - else { - i = findTime(timeStamp, *delayStruct); - assert(i < length); - time0 = (*delayStruct)[i].time; - value0 = (*delayStruct)[i].value; - // Was it the last value? - if (i+1 == length) { - if (i>0 && delayMax == delayTime) (*delayStruct).dequeue_n_first(i-1); - return value0; - } - time1 = (*delayStruct)[i+1].time; - value1 = (*delayStruct)[i+1].value; - if (i>0 && delayMax == delayTime) (*delayStruct).dequeue_n_first(i-1); + + if(delayTime < 0.0) + { + throw TerminateSimulationException(globalData->timeValue, + std::string("Negative delay requested.\n")); } - // was it an exact match? - if (time0 == timeStamp) { - //fprintf(stderr, "delayImpl: Exact match at %lf\n", currentTime); - return value0; + int length = ringBufferLength(delayStruct); + + if(length == 0) + { + // This occurs in the initialization phase + // fprintf(stderr, "delayImpl: Missing initial value, using argument value %g instead.\n", exprValue); + return exprValue; } - else if (time1 == timeStamp) { - return value1; + + // Returns: expr(time?delayTime) for time>time.start + delayTime and + // expr(time.start) for time <= time.start + delayTime. + // The arguments, i.e., expr, delayTime and delayMax, need to be subtypes of Real. + // DelayMax needs to be additionally a parameter expression. + // The following relation shall hold: 0 <= delayTime <= delayMax, + // otherwise an error occurs. If delayMax is not supplied in the argument list, + // delayTime need to be a parameter expression. See also Section 3.7.2.1. + // For non-scalar arguments the function is vectorized according to Section 10.6.12. + if(time <= tStart + delayTime) + { + double res = ((t_TimeAndValue*)getRingData(delayStruct, 0))->value; + // fprintf(stderr, "findTime: time <= tStart + delayTime: [%d] = %g\n",exprNumber,res); + return res; } - else { - //fprintf(stderr, "delayImpl: Linear interpolation of %lf between %lf and %lf\n", timeStamp, time0, time1); - - // linear interpolation - double timedif = time1 - time0; - double dt0 = time1 - timeStamp; - double dt1 = timeStamp - time0; - return (value0 * dt0 + value1 * dt1) / timedif; + else + { + // return expr(time-delayTime) + assert(delayTime >= 0.0); + double timeStamp = time - delayTime; + double time0, time1, value0, value1; + + int i; + + // find the row for the lower limit + if(timeStamp > ((t_TimeAndValue*)getRingData(delayStruct, length - 1))->time) + { + // delay between the last accepted time step and the current time + time0 = ((t_TimeAndValue*)getRingData(delayStruct, length - 1))->time; + value0 = ((t_TimeAndValue*)getRingData(delayStruct, length - 1))->value; + time1 = time; + value1 = exprValue; + } + else + { + i = findTime(timeStamp, *delayStruct); + assert(i < length); + time0 = ((t_TimeAndValue*)getRingData(delayStruct, i))->time; + value0 = ((t_TimeAndValue*)getRingData(delayStruct, i))->value; + + // Was it the last value? + if(i+1 == length) + { + if(i>0 && delayMax == delayTime) + dequeueNFirstRingDatas(delayStruct, i-1); + return value0; + } + time1 = ((t_TimeAndValue*)getRingData(delayStruct, i+1))->time; + value1 = ((t_TimeAndValue*)getRingData(delayStruct, i+1))->value; + if(i>0 && delayMax == delayTime) + dequeueNFirstRingDatas(delayStruct, i-1); + } + // was it an exact match? + if(time0 == timeStamp) + { + //fprintf(stderr, "delayImpl: Exact match at %lf\n", currentTime); + return value0; + } + else if(time1 == timeStamp) + { + return value1; + } + else + { + //fprintf(stderr, "delayImpl: Linear interpolation of %lf between %lf and %lf\n", timeStamp, time0, time1); + + // linear interpolation + double timedif = time1 - time0; + double dt0 = time1 - timeStamp; + double dt1 = timeStamp - time0; + return (value0 * dt0 + value1 * dt1) / timedif; + } } - } } diff --git a/c_runtime/Makefile.common b/c_runtime/Makefile.common index 832e91d2ad3..3b2cecf9eab 100644 --- a/c_runtime/Makefile.common +++ b/c_runtime/Makefile.common @@ -32,7 +32,7 @@ simulation_varinfo.o read_matlab4.o read_csv.o \ simulation_init.o simulation_input.o simulation_input_xml.o \ simulation_events.o linearize.o solver_main.o \ simulation_result_plt.o simulation_result_csv.o \ -simulation_result_mat.o simulation_delay.o tables.o options.o dgesv_aux.o \ +simulation_result_mat.o ringbuffer.o simulation_delay.o tables.o options.o dgesv_aux.o \ simulation_modelinfo.o \ $(META_OBJS) \ $(EXTRA_OBJS) \ diff --git a/c_runtime/ringbuffer.c b/c_runtime/ringbuffer.c new file mode 100644 index 00000000000..c76e18a21d1 --- /dev/null +++ b/c_runtime/ringbuffer.c @@ -0,0 +1,94 @@ +/* + * This file is part of OpenModelica. + * + * Copyright (c) 1998-2008, Linköpings University, + * Department of Computer and Information Science, + * SE-58183 Linköping, Sweden. + * + * All rights reserved. + * + * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THIS OSMC PUBLIC + * LICENSE (OSMC-PL). ANY USE, REPRODUCTION OR DISTRIBUTION OF + * THIS PROGRAM CONSTITUTES RECIPIENT'S ACCEPTANCE OF THE OSMC + * PUBLIC LICENSE. + * + * The OpenModelica software and the Open Source Modelica + * Consortium (OSMC) Public License (OSMC-PL) are obtained + * from Linköpings University, either from the above address, + * from the URL: http://www.ida.liu.se/projects/OpenModelica + * and in the OpenModelica distribution. + * + * This program is distributed WITHOUT ANY WARRANTY; without + * even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH + * IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS + * OF OSMC-PL. + * + * See the full OSMC Public License conditions for more details. + * + */ + +#include "ringbuffer.h" + +#include +#include +#include + +void allocRingBuffer(RINGBUFFER* rb, int sz, int item_size) +{ + rb->first_element = 0; + rb->num_element = 0; + rb->buf_size = sz > 0 ? sz : 1; + rb->item_size = item_size; + rb->buffer = calloc(rb->buf_size, rb->item_size); + assert(rb->buffer); +} + +void freeRingBuffer(RINGBUFFER* rb) +{ + free(rb->buffer); +} + +void* getRingData(RINGBUFFER* rb, int nIndex) +{ + assert(nIndex < rb->num_element); + assert(0 <= nIndex); + return ((char*)rb->buffer)+(((rb->first_element+nIndex)%rb->buf_size)*rb->item_size); +} + +void expandRingBuffer(RINGBUFFER* rb) +{ + int i; + void* temp = calloc(2*rb->buf_size, rb->item_size); + assert(temp); + + for(i=0; inum_element; i++) + memcpy(((char*)temp)+(i*rb->item_size), getRingData(rb, i), rb->item_size); + + free(rb->buffer); + rb->buffer = temp; + rb->buf_size *= 2; + rb->first_element = 0; +} + +void appendRingData(RINGBUFFER* rb, void* value) +{ + if(rb->buf_size < rb->num_element+1) + expandRingBuffer(rb); + + memcpy(((char*)rb->buffer)+(((rb->first_element+rb->num_element)%rb->buf_size)*rb->item_size), value, rb->item_size); + ++rb->num_element; +} + +void dequeueNFirstRingDatas(RINGBUFFER* rb, int n) +{ + assert(n <= rb->num_element); + assert(0 <= n); + rb->first_element = (rb->first_element+n)%rb->buf_size; + rb->num_element -= n; +} + +int ringBufferLength(RINGBUFFER* rb) +{ + return rb->num_element; +} diff --git a/c_runtime/ringbuffer.h b/c_runtime/ringbuffer.h index 9c3a3b4a1d4..86a1f005828 100644 --- a/c_runtime/ringbuffer.h +++ b/c_runtime/ringbuffer.h @@ -31,70 +31,37 @@ #ifndef RINGBUFFER_H #define RINGBUFFER_H -template - /* * This is an expanding ring buffer. * When it gets full, it doubles in size. * It's basically a queue which has get(ix) instead of get_first()/delete_first(). */ -class ringbuffer { -private: - - T* buffer; - int first_element; - int num_element; - int buf_size; - - T& fast_get(const int nIndex) { - return buffer[(first_element+nIndex)%buf_size]; - } - - void expand_buffer() { - T* nb = new T[buf_size*2]; - for (int i=0; i +#include "simulation_runtime.h" +#include "ringbuffer.h" #include "assert.h" #include "string.h" -#include "simulation_runtime.h" - +#include #include -#include "ringbuffer.h" double tStart = 0; typedef struct _TimeAndValue { - double time; - double value; + double time; + double value; } t_TimeAndValue; typedef struct _ExpressionDelayBuffer { - long currentIndex; - long maxExpressionBuffer; - t_TimeAndValue *expressionDelayBuffer; + long currentIndex; + long maxExpressionBuffer; + t_TimeAndValue *expressionDelayBuffer; } t_ExpressionDelayBuffer; // the delayStructure looks like a matrix (rows = expressionNumber+currentColumnIndex, columns={time, value}) -typedef ringbuffer t_buffer; -t_buffer **delayStructure; +RINGBUFFER **delayStructure; extern "C" { -extern const int numDelayExpressionIndex; + extern const int numDelayExpressionIndex; } + void initDelay(double startTime) { - // get the start time of the simulation: time.start. - tStart = startTime; - // allocate the memory for rows - delayStructure = new t_buffer*[numDelayExpressionIndex]; - for (int i=0; i= LOG_SOLVER) - { - fprintf(stdout, "initDelay called with startTime = %f\n", startTime); fflush(NULL); - } + int i; + + // get the start time of the simulation: time.start. + tStart = startTime; + + // allocate the memory for rows + delayStructure = new RINGBUFFER*[numDelayExpressionIndex]; + + for(i=0; i= LOG_SOLVER) + { + fprintf(stdout, "initDelay called with startTime = %f\n", startTime); + fflush(NULL); + } } void deinitDelay() { - delete[] delayStructure; + int i; + + for(i=0; i time) - end = i; - else - start = i; - } while (t != time && end > start + 1); - return start; + int start = 0; + int end = ringBufferLength(&delayStruct); + double t; + do + { + int i = (start + end) / 2; + t = ((t_TimeAndValue*)getRingData(&delayStruct, i))->time; + if(t > time) + end = i; + else + start = i; + }while (t != time && end > start + 1); + return start; } void storeDelayedExpression(int exprNumber, double exprValue) { - double time = globalData->timeValue; + t_TimeAndValue tpl; + double time = globalData->timeValue; + + if(exprNumber < 0) + { + // fprintf(stderr, "storeDelayedExpression: Invalid expression number %d.\n", exprNumber); + return; + } - if (exprNumber < 0) { - // fprintf(stderr, "storeDelayedExpression: Invalid expression number %d.\n", exprNumber); - return; - } + if(time < tStart) + { + // fprintf(stderr, "storeDelayedExpression: Time is smaller than starting time. Value ignored.\n"); + return; + } - if (time < tStart) { - // fprintf(stderr, "storeDelayedExpression: Time is smaller than starting time. Value ignored.\n"); - return; - } + // fprintf(stderr, "storeDelayed[%d] %g:%g\n", exprNumber, time, exprValue); - // fprintf(stderr, "storeDelayed[%d] %g:%g\n", exprNumber, time, exprValue); - - // Allocate more space for expressions - assert(exprNumber < numDelayExpressionIndex); + // Allocate more space for expressions + assert(exprNumber < numDelayExpressionIndex); - t_buffer* delayStruct = delayStructure[exprNumber]; - t_TimeAndValue tpl; - tpl.time = time; - tpl.value = exprValue; - delayStruct->append(tpl); + tpl.time = time; + tpl.value = exprValue; + appendRingData(delayStructure[exprNumber], &tpl); } double delayImpl(int exprNumber, double exprValue, double time, double delayTime, double delayMax /* Unused */) { - // fprintf(stderr, "delayImpl: exprNumber = %d, exprValue = %lf, time = %lf, delayTime = %lf\n", exprNumber, exprValue, time, delayTime); - - // Check for errors - assert(exprNumber >= 0); - assert(exprNumber < numDelayExpressionIndex); - - if (time <= tStart) { - // fprintf(stderr, "delayImpl: Entered at time < starting time: %g.\n", exprValue); - return exprValue; - } - - if (delayTime < 0.0) { - throw TerminateSimulationException(globalData->timeValue, - std::string("Negative delay requested.\n")); - } - - t_buffer* delayStruct = delayStructure[exprNumber]; - int length = delayStruct->length(); - - if (length == 0) { - // This occurs in the initialization phase - // fprintf(stderr, "delayImpl: Missing initial value, using argument value %g instead.\n", exprValue); - return exprValue; - } - - // Returns: expr(time?delayTime) for time>time.start + delayTime and - // expr(time.start) for time <= time.start + delayTime. - // The arguments, i.e., expr, delayTime and delayMax, need to be subtypes of Real. - // DelayMax needs to be additionally a parameter expression. - // The following relation shall hold: 0 <= delayTime <= delayMax, - // otherwise an error occurs. If delayMax is not supplied in the argument list, - // delayTime need to be a parameter expression. See also Section 3.7.2.1. - // For non-scalar arguments the function is vectorized according to Section 10.6.12. - - if (time <= tStart + delayTime) { - double res = (*delayStruct)[0].value; - // fprintf(stderr, "findTime: time <= tStart + delayTime: [%d] = %g\n",exprNumber,res); - return res; - } - else { - // return expr(time-delayTime) - assert(delayTime >= 0.0); - double timeStamp = time - delayTime; - double time0, time1, value0, value1; - - int i; + RINGBUFFER* delayStruct = delayStructure[exprNumber]; + + // fprintf(stderr, "delayImpl: exprNumber = %d, exprValue = %lf, time = %lf, delayTime = %lf\n", exprNumber, exprValue, time, delayTime); - // find the row for the lower limit - if (timeStamp > (*delayStruct)[length - 1].time) { - // delay between the last accepted time step and the current time - time0 = (*delayStruct)[length - 1].time; - value0 = (*delayStruct)[length - 1].value; - time1 = time; - value1 = exprValue; + // Check for errors + assert(0 <= exprNumber); + assert(exprNumber < numDelayExpressionIndex); + + if(time <= tStart) + { + // fprintf(stderr, "delayImpl: Entered at time < starting time: %g.\n", exprValue); + return exprValue; } - else { - i = findTime(timeStamp, *delayStruct); - assert(i < length); - time0 = (*delayStruct)[i].time; - value0 = (*delayStruct)[i].value; - // Was it the last value? - if (i+1 == length) { - if (i>0 && delayMax == delayTime) (*delayStruct).dequeue_n_first(i-1); - return value0; - } - time1 = (*delayStruct)[i+1].time; - value1 = (*delayStruct)[i+1].value; - if (i>0 && delayMax == delayTime) (*delayStruct).dequeue_n_first(i-1); + + if(delayTime < 0.0) + { + throw TerminateSimulationException(globalData->timeValue, + std::string("Negative delay requested.\n")); } - // was it an exact match? - if (time0 == timeStamp) { - //fprintf(stderr, "delayImpl: Exact match at %lf\n", currentTime); - return value0; + int length = ringBufferLength(delayStruct); + + if(length == 0) + { + // This occurs in the initialization phase + // fprintf(stderr, "delayImpl: Missing initial value, using argument value %g instead.\n", exprValue); + return exprValue; } - else if (time1 == timeStamp) { - return value1; + + // Returns: expr(time?delayTime) for time>time.start + delayTime and + // expr(time.start) for time <= time.start + delayTime. + // The arguments, i.e., expr, delayTime and delayMax, need to be subtypes of Real. + // DelayMax needs to be additionally a parameter expression. + // The following relation shall hold: 0 <= delayTime <= delayMax, + // otherwise an error occurs. If delayMax is not supplied in the argument list, + // delayTime need to be a parameter expression. See also Section 3.7.2.1. + // For non-scalar arguments the function is vectorized according to Section 10.6.12. + if(time <= tStart + delayTime) + { + double res = ((t_TimeAndValue*)getRingData(delayStruct, 0))->value; + // fprintf(stderr, "findTime: time <= tStart + delayTime: [%d] = %g\n",exprNumber,res); + return res; } - else { - //fprintf(stderr, "delayImpl: Linear interpolation of %lf between %lf and %lf\n", timeStamp, time0, time1); - - // linear interpolation - double timedif = time1 - time0; - double dt0 = time1 - timeStamp; - double dt1 = timeStamp - time0; - return (value0 * dt0 + value1 * dt1) / timedif; + else + { + // return expr(time-delayTime) + assert(delayTime >= 0.0); + double timeStamp = time - delayTime; + double time0, time1, value0, value1; + + int i; + + // find the row for the lower limit + if(timeStamp > ((t_TimeAndValue*)getRingData(delayStruct, length - 1))->time) + { + // delay between the last accepted time step and the current time + time0 = ((t_TimeAndValue*)getRingData(delayStruct, length - 1))->time; + value0 = ((t_TimeAndValue*)getRingData(delayStruct, length - 1))->value; + time1 = time; + value1 = exprValue; + } + else + { + i = findTime(timeStamp, *delayStruct); + assert(i < length); + time0 = ((t_TimeAndValue*)getRingData(delayStruct, i))->time; + value0 = ((t_TimeAndValue*)getRingData(delayStruct, i))->value; + + // Was it the last value? + if(i+1 == length) + { + if(i>0 && delayMax == delayTime) + dequeueNFirstRingDatas(delayStruct, i-1); + return value0; + } + time1 = ((t_TimeAndValue*)getRingData(delayStruct, i+1))->time; + value1 = ((t_TimeAndValue*)getRingData(delayStruct, i+1))->value; + if(i>0 && delayMax == delayTime) + dequeueNFirstRingDatas(delayStruct, i-1); + } + // was it an exact match? + if(time0 == timeStamp) + { + //fprintf(stderr, "delayImpl: Exact match at %lf\n", currentTime); + return value0; + } + else if(time1 == timeStamp) + { + return value1; + } + else + { + //fprintf(stderr, "delayImpl: Linear interpolation of %lf between %lf and %lf\n", timeStamp, time0, time1); + + // linear interpolation + double timedif = time1 - time0; + double dt0 = time1 - timeStamp; + double dt1 = timeStamp - time0; + return (value0 * dt0 + value1 * dt1) / timedif; + } } - } }