/* This is the version of the radau5 integrator without append and only for ODE */ #include #include "common.h" #include "corefunc.h" #include "stagen.h" #include "_linalg.h" #include "_util.h" #include "_integr2.h" /* Messages MACROS */ /* Return messges */ #define RetMessagem4 "Undefined user test function" #define RetMessagem3 "Matrix is repeatedly singular" #define RetMessagem2 "Step size become too small" #define RetMessagem11 "Round factor incosistent with inner constants" #define RetMessagem12 "Tolerances are too small or negative" #define RetMessagem13 "Number of newton iteration cannot be negative" #define RetMessagem14 "The number of indexed variables is inconsistent with the state space dimension" #define RetMessagem15 "Incosistent input for parameter m1 and/or m2" #define RetMessagem16 "Unacceptable value for step size control safety factor" #define RetMessagem17 "Unacceptable value for the jacobian computation cost" #define RetMessagem18 "Inconsistent newton stopping criterion and relative tolerance" #define RetMessagem19 "Inconsistent changing step size criterion parameter" #define RetMessagem110 "Inconsistent max and min step size ratio" #define RetMessagem111 "The bandwidth of mas cannot be smaller than thos of jac" #define RetMessagem112 "The hessenberg option can be applied only for ode with full jacobian" #define RetMessagem113 "The interval of integration cannot be negative" #define RetMessagem114 "The initial step cannot be negative" #define RetMessagem115 "The maximal step cannot be greater than the integraton interval" #define RetMessagem116 "UDF zero detection tolerance cannot be negative" #define RetMessagem117 "Maximum number of iteration for UDF zero detection cannot be negative" /* Error code for errors that do not allow the integrator to start */ #define NOSTARTERROR 4 /* UDF zeros messages */ #define UDFZeroMessage1 "Zero of user function number %i" #define UDFZeroMessage2 "Suspected zero of user function number %i" /* Integrator parameters MACROS */ #define _Tint contData->Tint /* Visible , it is set, out checked */ #define _h0 contData->h0 /* Visible , it is set, out checked */ #define _hmax contData->hmax /* Visible , it is set, out checked */ #define _hratiomin contData->hratiomin /* Invisible, to be set, inner checked */ #define _hratiomax contData->hratiomax /* Invisible, to be set, inner checked */ #define _chghratiomin contData->chghratiomin /* Invisible, to be set, inner checked */ #define _chghratiomax contData->chghratiomax /* Invisible, to be set, inner checked */ #define _hstrategy contData->hstrategy /* Invisible, to be set, out checked */ #define _hsafefact contData->hsafefact /* Invisible, to be set, inner checked */ #define _atol contData->atol /* Visible , it is set, inner checked */ #define _rtol contData->rtol /* Visible , it is set, inner checked */ #define _maxnewt contData->maxnewt /* Visible , it is set, inner checked */ #define _advnewt contData->advnewt /* Invisible, to be set, out checked */ #define _newtstop contData->newtstop /* Invisible, to be set, inner checked */ #define _ind1 contData->ind1 /* Invisible, to be set, inner checked */ #define _ind2 contData->ind2 /* Invisible, to be set, inner checked */ #define _ind3 contData->ind3 /* Invisible, to be set, inner checked */ #define _mljac contData->mljac /* Invisible, to be set, inner checked */ #define _mujac contData->mujac /* Invisible, to be set, inner checked */ #define _mlmas contData->mlmas /* Invisible, to be set, inner checked */ #define _mumas contData->mumas /* Invisible, to be set, inner checked */ #define _m1 contData->m1 /* Invisible, to be set, inner checked */ #define _m2 contData->m2 /* Invisible, to be set, inner checked */ #define _hess contData->hessenberg /* Invisible, to be set, inner checked */ #define _jaccost contData->jaccost /* Invisible, to be set, inner checked */ /* System class identifier */ #define _SystemCLASS contData->SystemCLASS /* Permanent */ /* UDF related MACROS*/ #define ActiveUTestFuns contData->ActiveUserTest #define nUTest genData->udFuncNum #define TestIter contData->IterTestFuns /* UDF maximum iteration for zero detection */ /* Visible , it is set, out checked */ #define _epsuf contData->epsuf /* UDF tolerance for zero detection */ /* Visible , it is set, out checked */ /* PostProcessor MACROS */ #define DefaultProcessor(p) ((ProcFuncPtr)genData->ProcFunc[0])(0,p,NULL,NULL) /* Messagge passing MACROS */ #define MsgProc(type,p,text) genData->ProcessPoint(type, p, text) /* Pointer type local variables, allow to access global data */ /*********************************************************** The pointer to Generator's functional parameters points to the following structure (names of functions must be listed in appropriate generator's description files under /genfunc section in proper order). ***********************************************************/ typedef struct { PrintVectorPtr PrintVector; CreateVectorPtr CreateVector; DestroyVectorPtr DestroyVector; GetVectorDimPtr GetVectorDim; CopyVectorPtr CopyVector; ScalProdPtr ScalProd; NormalizeVectorPtr NormalizeVector; AddScaledVectorPtr AddScaledVector; CopyVectorElementsPtr CopyVectorElements; } ContFuncPars, PNTR ContFuncParsPtr; /* Pointer to Data area passed by CONTENT */ Local(StagenDataPtr) genData=NULL; /* Pointer to a structure containing pointers to functional parameters */ Local(ContFuncParsPtr) contFunc=NULL; /* Pointer to a integrator's own data area */ Local(IntegrDataPtr) contData=NULL; /* Global copy of pointer to standard linea algebra functions */ Local(StdLinAlgDataPtr) staFunc=NULL; /* Global copy of pointer to utility functions */ Local(UtilitiesDataPtr) staUtil=NULL; /* More or less for all the subroutine */ Local(int) nphase, /* The state space dimension */ direction, /* The direction of the actual or previous integration */ ijac, /* Analitical or numerical jacobian 1/0 */ imas; /* ADE or ODE problem 1/0 */ Local(double) *X_Scratch=NULL; /* Phase vector for local purpouse */ Local(Vector) TX_Scratch=NULL, /* [T;X] vector needed for some functions */ TXnew_Scratch=NULL, /* [T;X] vector needed for some functions */ TXold_Scratch=NULL; /* [T;X] vector needed for some functions */ /* RHS jacket variables */ Local(Vector) F_Phase=NULL; /* The RHS values (just the ODE), the first component must be ignored it is (n+1)dim becouse in the integrators the state=[T;X]*/ /* Jac Jacket variables */ Local(Matrix) localJAC=NULL; /* Need to call the core Jacobian */ /* Mas Jacket variables */ Local(Matrix) localMAS=NULL; /* Need to call the core Mass */ /* UDF related variables, this variable are accessed by several routine */ Local(int) nmonitor, /* Number of monitored UDF */ *MonitorIndex=NULL, /* Indexes of monitored UDF */ ndetect, /* Number of detected UDF */ *DetectIndex=NULL; /* Indexes of detected UDF */ Local(Vector) TestValues, /* Vector of the UDF values */ *SuspiciousPoints; /* Point (T+X) where the detected UDFs annhilate, SupiciousPoints[i] is [T;X] where UDF[i]~=0 */ Local(int) *SortIndex, /* Index of zero UDF sorted respect to time, time(UDF[SortIndex[i]]~=0) <= time(UDF[SortIndex[i+1]]~=0) */ *UDFIter; /* UDFIter[i]=0 means Zero of UDF[i] detected within the prefixed number of dycotomy iteration UDFIter[i]=0 means that the iteration has been finishe before reaching the tollerance */ char MsgTxt[80]; /* Message string, to compose the out messages */ /* radau5.c working memory */ Local(double) *work=NULL; /* Real numbers work space */ Local(long int) *iwork=NULL; /* Integer numbers work space */ /* EntryPoint mode*/ typedef EntryPtr(Int2,ProcFuncPtr,(Int2 Ind, Vector xx, Vector vv, CharPtr *msg)); /* Prototype */ Local (void) RHS_Jacket(long int, double, double *, double *); Local (void) JAC_Jacket(long int, double, double *, double *, long int); Local (int) EveryStepFunction(long int, double, double, double *, double *); Local (void) MAS_Jacket(long int, double *, long int); Local (void) DeallocateMemory(); double contr5rad(long int, double x, double *); /****************************************************************************** LOOK AT THE END FOR LOCAL FUNCTIONS ******************************************************************************/ /******************* BEGIN OF Entry Point radau5_OrbitCont *******************/ /****************************************************************************** Entry point of the Radau5 integrator It check the calling parameters and at the end simply call the local function radau5. ******************************************************************************/ Entry(Int2) radau5_OrbitCont(StagenDataPtr stagen) { /* Local variables */ int returnmode, /* to check the return mode of the integrator */ i,j, /* some index is always usefull */ ndimp1; /* nphase+1, it is used often */ long int ipar[20]; /* Advanced integer parameters */ double StartTime,EndTime, /* Start and end time of integration */ ReachedTime, /* Where the integration is arrived */ NextEndTime, /* Up to where go if the computation is extended */ rpar[20]; /* Advanced real parameters */ /* First of all check if the call is right, stagen must be not NULL */ if (!stagen) { DeallocateMemory(); return((Int2)(1)); } /* Local (and more handable) copy of global pointers */ genData=stagen; contData=(IntegrDataPtr)stagen->GenPtr; contFunc=(ContFuncParsPtr)stagen->GenFunc; staFunc=(StdLinAlgDataPtr)stagen->StdLinAlg; staUtil=(UtilitiesDataPtr)stagen->Utilities; /* Detect the dimesion of the system */ nphase=staFunc->GetVectorDim(contData->X0); ndimp1=nphase; nphase--; /************************************************************************* Check the parameter that are not checked inside the integrator, namely: _epsuf, TestIter, Tint, h0, hmax, hstrategy, advnewt the last three are set in this routine so they are right by definition *************************************************************************/ returnmode=0; if (_Tint<=((double)(0))) returnmode=-113; else if (_h0<=((double)(0))) returnmode=-114; else if (_hmax>_Tint) returnmode=-115; else if (_epsuf<=((double)(0))) returnmode=-116; else if (TestIter<=0) returnmode=-117; /* Getting how many UDF has swoitched monitor or detect */ for (i=0,nmonitor=0; iCreateVector(ndimp1); /* (2) */ TXnew_Scratch=staFunc->CreateVector(ndimp1); /* (3) */ TXold_Scratch=staFunc->CreateVector(ndimp1); /* (4) */ /*For the RHS Jacket*/ F_Phase=staFunc->CreateVector(ndimp1); /* (5) */ /*For the JAC Jacket*/ /* it has been moved inside the jacobian check */ /* (6.a) */ /*For the MAS Jacket*/ /* it has been moved inside the DAE/ODE check */ /* (6.b) */ /*For monitored and detected UDF*/ if (nmonitor) MonitorIndex=MemNew(sizeof(int)*nmonitor); /* (7) */ if (ndetect) { DetectIndex=MemNew(sizeof(int)*ndetect); /* (8) */ SuspiciousPoints=MemNew(ndetect*sizeof(Vector *)); /* (9) */ for(i=0;iCreateVector(ndimp1); /* (10) */ SortIndex=MemNew(sizeof(int)*ndetect); /* (11) */ UDFIter=MemNew(sizeof(int)*ndetect); /* (12) */ } TestValues=staFunc->CreateVector(nUTest); /* (13) */ /* I the outer checked parameter are ok */ if (!returnmode) { /* Build up the monitored and detected UDF indexes */ if (nmonitor) for (i=0,j=0; iCreateMatrix(nphase,nphase); /* (6.b) */ /* Determine numerical or routine for the jacobian, if localJAC is needed it is allocated */ ijac=(genData->Der1)?1:0; if (ijac) localJAC=staFunc->CreateMatrix(nphase,nphase); /* (6.a) */ /* Setting the values of invisible parameter at standard values */ /**************************************************************************** This part of the code is needed just too see the corrispondance between the real name of parameters and their position in rpar and ipar. For the meaning of the various parameter take a look to the file radau5.c ****************************************************************************/ ipar[0]=_hess=0; /* No Hessenberg decomposition */ ipar[3]=_advnewt=0; /* Advanced initial guess is used */ ipar[7]=_hstrategy=0; /* Advanced step size strategy */ ipar[8]=_m1=0; /* Not a second order system */ ipar[9]=_m2=0; /* Not a second order system */ rpar[1]=_hsafefact=.0; rpar[2]=_jaccost=.0; rpar[3]=_newtstop=.0; rpar[4]=_chghratiomin=.0; rpar[5]=_chghratiomax=.0; rpar[7]=_hratiomin=.0; rpar[8]=_hratiomax=.0; /* Setting the visible advanced parameters */ ipar[2]=_maxnewt; ipar[4]=_ind1; ipar[5]=_ind2; ipar[6]=_ind3; rpar[6]=_hmax; /* THIS HAS NOT TO BE TOUCHED */ rpar[0]=((double)(5e-19)); /* Internal maximal resolution of radau5 */ /* Determining if is a new computation or an extent of the previous */ /********************************************************************* After the if or else branch StartTime EndTime and the initial state X_Scratch are set *********************************************************************/ /********************************************************************* WARNING: the follow assume that extend is allowed if and only if nothing has been changed since the last computation (DO NOT TOUCH THE PARAMETERS ... AND SO ON) *********************************************************************/ /*********************************************************************** Internal structure saved at genData->ContDataGenPtr [0] time of the last point of the previous integration (t) [1 ... nphase] state of the last point of the previous integration (X) [nphase+1] the pre-stored EndTime (Treached+Tint) [nphase+2] direction of previous integration (0=bw/1=fw) ***********************************************************************/ if (genData->ContDataGenLen) /* Extends the previous integration */ { /* Get the reached time of integration */ StartTime=(genData->ContDataGenPtr)[0]; /* Get the reached state */ MemCopy(X_Scratch,&(genData->ContDataGenPtr)[1],nphase*sizeof(VAR)); /* Get the pre-stored EndTime */ EndTime=genData->ContDataGenPtr[ndimp1]; /* Get the previous direction */ direction=(int)(genData->ContDataGenPtr[ndimp1+1]); /* Free the Data for extending, they will be allocated again at the end of the integration */ genData->ContDataGenPtr=MemFree(genData->ContDataGenPtr); } else /* Start a new one integration */ { /* Get the initial time */ StartTime=((FloatHiPtr)(contData->X0))[0]; /* Get the dirction */ direction=genData->Forward; /* Use the direction to compute the end time */ EndTime=direction ? StartTime+_Tint : StartTime-_Tint; /* Get the initial state */ MemCopy(X_Scratch,&((FloatHiPtr)(contData->X0))[1],nphase*sizeof(VAR)); /*********************************************** WARNING: who free the genData->ContDataGenPtr if this branch is chosen, this is equal to the other two intgrators but could be a bug ***********************************************/ } /* Now the integrator can be called */ /********************************************************** For a description of this function take a look at the radau5.c file **********************************************************/ returnmode=Radau5Integrator( nphase, RHS_Jacket, StartTime, EndTime, &ReachedTime, X_Scratch, _h0, _rtol, _atol, JAC_Jacket, ijac, _mljac, _mujac, MAS_Jacket, imas, _mlmas, _mumas, rpar,ipar, EveryStepFunction ); /* Compose the TX vector need for the post processing */ TX_Scratch[0]=ReachedTime; MemCopy(&TX_Scratch[1],X_Scratch,sizeof(VAR)*nphase); } /* Handle correctly the returned value */ switch(returnmode) { case -4: /* Computation of an UDF failed */ MsgProc(0,TX_Scratch,RetMessagem4); returnmode=2; /* Call to the default postprocessor */ DefaultProcessor(TX_Scratch); break; case -3: /* Matrix is repeatedly singular */ MsgProc(0,TX_Scratch,RetMessagem3); returnmode=1; /* Call to the default postprocessor */ DefaultProcessor(TX_Scratch); break; case -2: /* Step size become too small */ MsgProc(0,TX_Scratch,RetMessagem2); returnmode=3; /* Call to the default postprocessor */ DefaultProcessor(TX_Scratch); break; case -11: /* round factor incosistent with inner constants */ staUtil->ErrMessage(RetMessagem11); returnmode=NOSTARTERROR; break; case -12: /* tolerances are too small or negative */ staUtil->ErrMessage(RetMessagem12); returnmode=NOSTARTERROR; break; case -13: /* number of newton iteration cannot be negative */ staUtil->ErrMessage(RetMessagem13); returnmode=NOSTARTERROR; break; case -14: /* the number of indexed variables is inconsistent with the dimension n */ staUtil->ErrMessage(RetMessagem14); returnmode=NOSTARTERROR; break; case -15: /* incosistent input for parameter m1 and/or m2 */ staUtil->ErrMessage(RetMessagem15); returnmode=NOSTARTERROR; break; case -16: /* unacceptable value for step size control safety factor */ staUtil->ErrMessage(RetMessagem16); returnmode=NOSTARTERROR; break; case -17: /* unacceptable value for the jacobian computation cost */ staUtil->ErrMessage(RetMessagem17); returnmode=NOSTARTERROR; break; case -18: /* inconsistent newton stopping criterion and relative tolerance */ staUtil->ErrMessage(RetMessagem18); returnmode=NOSTARTERROR; break; case -19: /* inconsistent changing step size criterion parameter */ staUtil->ErrMessage(RetMessagem19); returnmode=NOSTARTERROR; break; case -110: /* inconsistent max and min step size ratio */ staUtil->ErrMessage(RetMessagem110); returnmode=NOSTARTERROR; break; case -111: /* the bandwidth of mas cannot be smaller than thos of jac */ staUtil->ErrMessage(RetMessagem111); returnmode=NOSTARTERROR; break; case -112: /* the hessenberg option can be applied only for ode with full jacobian */ staUtil->ErrMessage(RetMessagem112); returnmode=NOSTARTERROR; break; case -113: /* The interval of integration cannot be negative */ staUtil->ErrMessage(RetMessagem113); returnmode=NOSTARTERROR; break; case -114: /* The initial step cannot be negative */ staUtil->ErrMessage(RetMessagem114); returnmode=NOSTARTERROR; break; case -115: /* The maximal step cannot be greater than the integraton interval */ staUtil->ErrMessage(RetMessagem115); returnmode=NOSTARTERROR; break; case -116: /* UDF zero detection tolerance is negative */ staUtil->ErrMessage(RetMessagem116); returnmode=NOSTARTERROR; break; case -117: /* Maximum number of iteration for UDF zero detection is negative */ staUtil->ErrMessage(RetMessagem117); returnmode=NOSTARTERROR; break; case 0: /* The computation has been ended normally */ NextEndTime=direction ? EndTime+_Tint : EndTime-_Tint; returnmode=0; /*this is stupid but is just for simmetry */ /* Call to the default postprocessor */ DefaultProcessor(TX_Scratch); break; case 1: /* The computation has been stopped normally */ NextEndTime=EndTime; returnmode=0; /* Call to the default postprocessor */ DefaultProcessor(TX_Scratch); break; } /* This part has to be done only if the integrator has done at least one step */ if (!returnmode) { /* Storing the total or partial result of the computation for the eventual extend */ /*********************************************************************** Internal structure saved at genData->ContDataGenPtr [0] time of the last point of the previous integration (t) [1 ... nphase] state of the last point of the previous integration (X) [nphase+1] the pre-stored EndTime (Treached+Tint) [nphase+2] direction of previous integration (0=bw/1=fw) ***********************************************************************/ genData->ContDataGenLen=(ndimp1+2)*sizeof(VAR); genData->ContDataGenPtr=MemNew(genData->ContDataGenLen); MemCopy(genData->ContDataGenPtr,TX_Scratch,ndimp1*sizeof(VAR)); genData->ContDataGenPtr[ndimp1]=NextEndTime; genData->ContDataGenPtr[ndimp1+1]=(VAR)(direction); } /* Deallocate the memory */ DeallocateMemory(); /* Return with the correct return code */ return((Int2)(returnmode)); } /********************* END OF Entry Point radau5_OrbitCont *******************/ /********************** BEGIN OF DeallocateMemory ****************************/ /****************************************************************************** ******************************************************************************/ Local (void) DeallocateMemory(void) { int i; X_Scratch=MemFree(X_Scratch); /* (1) */ TX_Scratch=staFunc->DestroyVector(TX_Scratch); /* (2) */ TXnew_Scratch=staFunc->DestroyVector(TXnew_Scratch); /* (3) */ TXold_Scratch=staFunc->DestroyVector(TXold_Scratch); /* (4) */ F_Phase=staFunc->DestroyVector(F_Phase); /* (5) */ if (ijac) localJAC=staFunc->DestroyMatrix(localJAC); /* (6.a) */ if (imas) localMAS=staFunc->DestroyMatrix(localMAS); /* (6.b) */ if (nmonitor) MonitorIndex=MemFree(MonitorIndex); /* (7) */ if (ndetect) { for(i=0;iDestroyVector(SuspiciousPoints[i]);/* (9) */ SuspiciousPoints=MemFree(SuspiciousPoints); /* (8) */ DetectIndex=MemFree(DetectIndex); /* (10) */ SortIndex=MemFree(SortIndex); /* (11) */ UDFIter=MemFree(UDFIter); /* (12) */ } TestValues=staFunc->DestroyVector(TestValues); /* (13) */ iwork=MemFree(iwork); work=MemFree(work); } /************************* END OF DeallocateMemory ***************************/ /*************************** BEGIN OF RHS_Jacket *** *************************/ /****************************************************************************** Local (void) RHS_Jacket(long int n,double t,double *X,double *DXDT) It is the local jacket for the RHS function that allows radau5 to evaluate the RHS. WARNING: there is no check on the failure of the computation of RHS. This is because the integrator do not provide a way to check that. This has to be modified in such a way that the integrator go out if the RHS fail in the evaluation. Anyway, at the present shape the CONTENT core function that evaluate the RHS return always 0 so the check is unusefull. Analogous warning is valid for the jacobian too. ******************************************************************************/ Local (void) RHS_Jacket(long int n,double t,double *X,double *DXDT) { /* Compose the TX vector */ TX_Scratch[0]=t; MemCopy(&TX_Scratch[1],X,sizeof(VAR)*nphase); /* Evaluating the RHS */ (genData->DefFunc[0])(TX_Scratch,F_Phase); /* Copy the RHS in the norval array DXDT */ MemCopy(DXDT,&F_Phase[1],nphase*sizeof(VAR)); } /******************************** END OF RHSJacket ***************************/ /*************************** BEGIN OF JacobianJacket *************************/ /****************************************************************************** Local (void) JAC_Jacket(long int n,double t,double *X,double *DFDX, long int rowDFDX) It is the local jacket for the JAC function that allows radau5 to evaluate the JAC itself. WARNING: there is no check on the failure of the computation of JAC. This is because the integrator do not provide a way to check that. This has to be modified in such a way that the integrator go out if the JAC fail in the evaluation. Anyway, at the present shape the CONTENT core function that evaluate the JAC return always 0 so the check is unusefull. Analogous warning is valid for the RHS too. For the moment rowDFDX is not used because the jacobian is always assumed to be full, it can be extended in the future. WARNING: In order to let it working, __o.c and __o.con has been changed because it was not available the function DefFunc[1] for the orbit continuation. ******************************************************************************/ Local (void) JAC_Jacket(long int n,double t,double *X,double *DFDX, long int rowDFDX) { int i,j,ll,ul; /* Some local index */ /* Compose the TX vector */ TX_Scratch[0]=t; MemCopy(&TX_Scratch[1],X,sizeof(VAR)*nphase); /* Evaluating the Jacobian */ (genData->DefFunc[1])(TX_Scratch,localJAC); /* Copy the Jacobian from the matrix localJAC to DFDX */ if(!(n-rowDFDX)) /*Full Jacobian*/ for (i=0;i0)?(i-_mljac):0; ul=((i+_mujac)DefFunc[2])(localMAS[0]); /* Copy the Mass Matrix from the matrix localMAS to M */ if(!(n-rowM)) /*Full MasMatrix*/ for (i=0;i0)?(i-_mlmas):0; ul=((i+_mumas)udFunc)(TXnew_Scratch,MonitorIndex[i],&TestValues[MonitorIndex[i]]); for (i=0; (iudFunc)(TXnew_Scratch,DetectIndex[i],&TestValues[DetectIndex[i]]); } else /* This is not the first point */ { /* Compose the TX vector in the new point */ TXnew_Scratch[0]=tn; MemCopy(&TXnew_Scratch[1],X,sizeof(VAR)*nphase); /* Compose the TX vector in the old point */ TXold_Scratch[0]=tnm1; for (tmp=1;tmp<=n;tmp++) TXold_Scratch[tmp]=contr5rad(tmp,TXold_Scratch[0],CV); /* Check and detect eventually zeros of UDF in the interval tnm1-tn */ BadTestFunctions=0; /* Up to now no UDF evaluation */ UDFZero=0; /* No zeros detected yet */ /* This is for the detect UDF */ for (i=0;(iudFunc)(TXnew_Scratch,DetectIndex[i],&TstVal2); /* Value of the test function in the old point, previously saved */ TstVal1=TestValues[DetectIndex[i]]; /* Save the current value */ TestValues[DetectIndex[i]]=TstVal2; /* Check if in the first extreme is zero */ if (TstVal1==((double)(0))) { UDFZero=1; /* At least one zeros detected */ SortIndex[i]=i; /* Later it will be sorted */ UDFIter[i]=0; /* Detection inside the iteration, in such case none */ SuspiciousPoints[SortIndex[i]]=contFunc->CopyVector(TXold_Scratch,SuspiciousPoints[SortIndex[i]]); } /******************************************************************** The check if in the second extreme is zero is not necessary because it will be done at the next step ********************************************************************/ else if (TstVal1*TstVal2 < 0.0) /* A zero must exist between Told and Tnew */ { UDFZero=1; /* At least one zeros detected */ SortIndex[i]=i; /* Later it will be sorted */ /* Dicotomy on the interval Told Tnew */ tau1=tnm1; tau2=tn; outsidetol=1; /* Look for the zero by dicotomy with at maximum TestIter Iteration */ for (iter=0; (iterudFunc)(TX_Scratch,DetectIndex[i],&TstVal); /* Dicotomy step */ if(!BadTestFunctions) { if (fabs(TstVal)<_epsuf) /* The zero has been found */ outsidetol=0; else if (TstVal*TstVal1 < 0.0) /* New interval is tau1-tau */ { TstVal2=TstVal; tau2=tau; } else /* New interval is tau-tau2 */ { TstVal1=TstVal; tau1=tau; } } } /****************************************************************************************** The version in RungeKutta and Euler there is no a paricular action for iter>TestIter, this does not really make sense. In this implemntation the zeros are subdivided in suspected zero (should be there a zero but has not be found with TestIter iteration) zero (is there a zero and has been found within TestIter iteration) ******************************************************************************************/ if (iter>=TestIter) UDFIter[i]=1; /* Detection not inside the iteration */ else UDFIter[i]=0; /* Detection inside the iteration */ /* Save the zero point */ SuspiciousPoints[SortIndex[i]]=contFunc->CopyVector(TX_Scratch,SuspiciousPoints[SortIndex[i]]); } } /* Sorting of zeros respect to time X[0] */ if (UDFZero&&(!BadTestFunctions)) { for (i=0;i=0)&&(SortIndex[j]>=0)) /* The i-th and j-th detected test functions are zero */ if (direction ? (SuspiciousPoints[SortIndex[i]][0]>SuspiciousPoints[SortIndex[j]][0]) : (SuspiciousPoints[SortIndex[i]][0]=0) if (!UDFIter[i]) /* Detected zero */ { sprintf(MsgTxt,UDFZeroMessage1,DetectIndex[SortIndex[i]]+1); if (MsgProc(USERPOINT|(DetectIndex[SortIndex[i]]),SuspiciousPoints[SortIndex[i]],MsgTxt)) return(1); } else /* Suspected zero */ { sprintf(MsgTxt,UDFZeroMessage2,DetectIndex[SortIndex[i]]+1); if (MsgProc(USERPOINT|(DetectIndex[SortIndex[i]]),SuspiciousPoints[SortIndex[i]],MsgTxt)) return(1); } } /* Recompute values of all test function at the next regular point, alias Tnew */ for (i=0; (iudFunc)(TXnew_Scratch,MonitorIndex[i],&TestValues[MonitorIndex[i]]); for (i=0; (iudFunc)(TXnew_Scratch,DetectIndex[i],&TestValues[DetectIndex[i]]); } /* Storing and Stop-check Part */ if (BadTestFunctions) return (2); else return(MsgProc(0,TXnew_Scratch,NULL)?1:0); } /************************** END OF EveryStepFunction *************************/ /* Inclusion of the radau library file */ #include"_radau5.c"