/[ascend]/trunk/tcltk/generic/interface/Integrators.c
ViewVC logotype

Diff of /trunk/tcltk/generic/interface/Integrators.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 571 by johnpye, Tue May 9 00:14:59 2006 UTC revision 670 by johnpye, Wed Jun 21 07:57:19 2006 UTC
# Line 1  Line 1 
1  /*  /*  ASCEND modelling environment
2   *  Integrators.c      Copyright 1997, Carnegie Mellon University
3   *  by Kirk Abbott and Ben Allan      Copyright (C) 2006 Carnegie Mellon University
4   *  Created: 1/94  
5   *  Version: $Revision: 1.32 $      This program is free software; you can redistribute it and/or modify
6   *  Version control file: $RCSfile: Integrators.c,v $      it under the terms of the GNU General Public License as published by
7   *  Date last modified: $Date: 2003/08/23 18:43:06 $      the Free Software Foundation; either version 2, or (at your option)
8   *  Last modified by: $Author: ballan $      any later version.
9   *  
10   *  This file is part of the ASCEND Tcl/Tk interface      This program is distributed in the hope that it will be useful,
11   *      but WITHOUT ANY WARRANTY; without even the implied warranty of
12   *  Copyright 1997, Carnegie Mellon University      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13   *      GNU General Public License for more details.
14   *  The ASCEND Tcl/Tk interface is free software; you can redistribute  
15   *  it and/or modify it under the terms of the GNU General Public License as      You should have received a copy of the GNU General Public License
16   *  published by the Free Software Foundation; either version 2 of the      along with this program; if not, write to the Free Software
17   *  License, or (at your option) any later version.      Foundation, Inc., 59 Temple Place - Suite 330,
18   *      Boston, MA 02111-1307, USA.
19   *  The ASCEND Tcl/Tk interface is distributed in hope that it will be  *//**
20   *  useful, but WITHOUT ANY WARRANTY; without even the implied warranty of      @file
21   *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU      Tcl/Tk interface functions for the Integration feature
22   *  General Public License for more details.  *//*
23   *      by Kirk Abbott and Ben Allan
24   *  You should have received a copy of the GNU General Public License      Created: 1/94
25   *  along with the program; if not, write to the Free Software Foundation,      Version: $Revision: 1.32 $
26   *  Inc., 675 Mass Ave, Cambridge, MA 02139 USA.  Check the file named      Version control file: $RCSfile: Integrators.c,v $
27   *  COPYING.  COPYING is found in ../compiler.      Date last modified: $Date: 2003/08/23 18:43:06 $
28   */      Last modified by: $Author: ballan $
29    */
30  /*  
31   * This module defines the general integration auxillaries for Ascend.  #include <tcl.h>
32   */  
33    #include <utilities/ascConfig.h>
34  #include <time.h>  #include <solver/integrator.h>
35  #include <tcl.h>  #include <solver/lsode.h>
36  #include <utilities/ascConfig.h>  #include <solver/samplelist.h>
37  #include <utilities/ascMalloc.h>  
38  #include <utilities/readln.h>  #include "HelpProc.h"
39  #include <general/list.h>  #include "Integrators.h"
40  #include <general/dstring.h>  #include "BrowserQuery.h"
41  #include <compiler/compiler.h>  #include "Qlfdid.h"
42  #include <compiler/instance_enum.h>  #include "UnitsProc.h"
43  #include <compiler/fractions.h>  #include "BrowserProc.h"
44  #include <compiler/dimen.h>  #include "HelpProc.h"
45  #include <compiler/units.h>  #include "SolverGlobals.h"
46  #include <compiler/module.h>  
47  #include <compiler/library.h>  #ifndef lint
48  #include <compiler/types.h>  static CONST char IntegratorsID[] = "$Id: Integrators.c,v 1.32 2003/08/23 18:43:06 ballan Exp $";
49  #include <compiler/child.h>  #endif
50  #include <compiler/type_desc.h>  
51  #include <compiler/atomvalue.h>  #define SNULL (char *)NULL
52  #include <compiler/instance_name.h>  
53  #include <compiler/instquery.h>  static SampleList l_samplelist;
54  #include <compiler/parentchild.h>  
55  #include <compiler/symtab.h>  /*-------------------------------------------------------
56  #include <compiler/instance_io.h>    HANDLING OF OUTPUT FILES
57  #include <solver/slv_types.h>  */
58  #include <solver/mtx.h>  
59  #include <solver/var.h>  /* vars relating to output files */
60  /*  static FILE *l_obs_file = NULL;
61   * note: the analytic jacobian routines (state matrix) depend on the  static char *l_obs_filename = NULL;
62   * assumption that struct var_variable *<--> struct Instance *.  
63   */  static FILE *l_y_file = NULL;
64  #include <solver/rel.h>  static char *l_y_filename = NULL;
65  #include <solver/discrete.h>  
66  #include <solver/conditional.h>  static int l_print_option = 1; /* default si */
67  #include <solver/logrel.h>  static int l_print_fmt = 0; /* default variable */
68  #include <solver/bnd.h>  
69  #include <solver/slv_common.h>  FILE *Asc_IntegOpenYFile(void)
70  #include <solver/linsol.h>  {
71  #include <solver/linsolqr.h>    if (l_y_filename==NULL) {
72  #include <solver/slv_client.h>      return NULL;
73  #include "HelpProc.h"    }
74  #include "Integrators.h"    l_y_file = fopen(l_y_filename,"a+");
75  #include "BrowserQuery.h"  
76  #include "Qlfdid.h"    if (l_y_file==NULL) {
77  #include "UnitsProc.h"      FPRINTF(ASCERR,
78  #include "BrowserProc.h"        "WARNING: (integrate) Unable to open\n\t%s\nfor state output log.\n",
79  #include "HelpProc.h"        l_y_filename);
80  #include "SolverGlobals.h"    } else {
81  #include "Lsode.h"      time_t t;
82        t = time((time_t *)NULL);
83  #ifndef lint      FPRINTF(l_y_file,"DATASET %s", asctime(localtime(&t)));
84  static CONST char IntegratorsID[] = "$Id: Integrators.c,v 1.32 2003/08/23 18:43:06 ballan Exp $";      FFLUSH(l_y_file);
85  #endif    }
86      return l_y_file;
87    }
88  #define SNULL (char *)NULL  
89    FILE *Asc_IntegOpenObsFile(void)
90  static symchar *g_symbols[3];  {
91      if (l_obs_filename==NULL) {
92  /* The following names are of solver_var children or attributes      return NULL;
93   * we support (at least temporarily) to determine who is a state and    }
94   * who matching derivative.    l_obs_file = fopen(l_obs_filename,"a+");
95   * These should be supported directly in a future solveratominst.    if (l_obs_file==NULL) {
96   */      FPRINTF(ASCERR,
97  #define STATEFLAG g_symbols[0]        "WARNING: (integrate) Unable to open\n\t%s\nfor observation log.\n",
98  /* Integer child. 0= algebraic, 1 = state, 2 = derivative, 3 = 2nd deriv etc */        l_obs_filename);
99  #define STATEINDEX g_symbols[1]    } else {
100  /* Integer child. all variables with the same STATEINDEX value are taken to      time_t t;
101   * be derivatives of the same state variable. We really need a compiler      t = time((time_t *)NULL);
102   * that maintains this info by backpointers, but oh well until then.      FPRINTF(l_obs_file,"DATASET %s", asctime(localtime(&t)));
103   */      FFLUSH(l_obs_file);
104  #define OBSINDEX g_symbols[2]    }
105  /* Integer child. All variables with OBSINDEX !=0 will be recorded in    return l_obs_file;
106   * the blsode output file. Tis someone else's job to grok this output.  }
107   */  
108    FILE *Asc_IntegGetYFile(void)
109  #define INTDEBUG 0  {
110      return l_y_file;
111  struct Integ_var_t {  }
112    long index;  
113    long type;  FILE *Asc_IntegGetObsFile(void)
114    struct var_variable *i;  {
115  };    return l_obs_file;
116  /* temp catcher of dynamic variable and observation variable data */  }
117    
118  struct Integ_samples_t {  void Asc_IntegReleaseYFile(void)
119    long ns;  {
120    double *x;    l_y_file = NULL;
121    dim_type *d;  }
122  };  void Asc_IntegReleaseObsFile(void)
123  /* an array of sample 'times' for reference during integration */  {
124      l_obs_file = NULL;
125  /* global to the world */  }
126  struct Integ_DiffData g_intg_diff = {0L,NULL,NULL,NULL,NULL,NULL,NULL};  
127    /********************************************************************/
128  /*Define items for integration interface. they are global to this file.  /* string column layout: 26 char columns  if l_print_fmt else variable */
129  this globalness is bad. */  /* numeric format: leading space plus characters from Unit*Value */
130  static slv_system_t g_intgsys_cur = NULL;  /* index format: left padded long int */
131      /* derivative system children */  /* headline format space plus - s */
132  static struct var_variable *g_intginst_d_x = NULL;  #define BCOLSFMT (l_print_fmt ? "%-26s" : "\t%s")
133      /* derivative system sizes */  #define BCOLNFMT (l_print_fmt ? " %-25s" : "\t%s")
134  static long g_intginst_d_n_eq, g_intginst_d_n_obs;  #define BCOLIFMT (l_print_fmt ? " %25ld" : "\t%ld")
135      /* derivative system array var lists. Not null terminated */  #define BCOLHFMT (l_print_fmt ? " -------------------------" : "\t---")
136  static struct var_variable **g_intginst_d_y_vars = NULL;  
137  static struct var_variable **g_intginst_d_dydx_vars = NULL;  void Asc_IntegPrintYHeader(FILE *fp, IntegratorSystem *blsys)
138  static struct var_variable **g_intginst_d_obs_vars = NULL;  {
139  /* sampling vector for feeding smarter integrators */    long i,len, *yip;
140  static struct Integ_samples_t g_intg_samples = {0L, NULL, NULL};    char *name;
141      struct Instance *in;
142      /* Integ_system_build locally global var */    int si;
143  static struct Integ_system_t *l_isys = NULL;    if (fp==NULL) {
144  static long l_nstates = 0;      return;
145  static long l_nderivs = 0;    }
146      if (blsys==NULL) {
147  /* Macros for d. arrays.      FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintYHeader: called w/o data\n");
148     D_Y_VAR(i) returns the atom that corresponds to the FORTRAN y(i).      return;
149     D_DYDX_VAR(i) returns the atom that corresponds to the FORTRAN ydot(i).    }
150     D_OBS_VAR(i) returns the atom that corresponds to the ASCEND d.obs[i].    if (blsys->n_y == 0) {
151  */      return;
152  #define D_Y_VAR(i) g_intginst_d_y_vars[(i)-1]    }
153  #define D_DYDX_VAR(i) g_intginst_d_dydx_vars[(i)-1]    if (blsys->y == NULL) {
154  #define D_OBS_VAR(i) g_intginst_d_obs_vars[(i)-1]      FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintYHeader: called w/NULL data\n");
155        return;
156  /********************************************************************/    }
157      len = blsys->n_y;
158  int Asc_IntegCheckStatus(slv_status_t status) {    yip = blsys->y_id;
159    if (status.converged) {    si = l_print_option;
160      return 1;    /* output indep var name */
161    }    /* output dep var names */
162    if (status.diverged) {    FPRINTF(fp,"States: (user index) (name) (units)\n");
163      FPRINTF(stderr, "The derivative system did not converge.\n");    in = var_instance(blsys->x);
164      FPRINTF(stderr, "Integration will be terminated ");    FPRINTF(fp,"{indvar}");  /* indep id name */
165      FPRINTF(stderr, "at the end of the current step.\n");    name = WriteInstanceNameString(in,g_solvinst_cur);
166      return 0;    FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));
167    }    ascfree(name);
168    if (status.inconsistent) {    for (i=0; i< len; i++) {
169      FPRINTF(stderr, "A numerical inconsistency was discovered ");      in = var_instance(blsys->y[i]);
170      FPRINTF(stderr, "while calculating derivatives.");      FPRINTF(fp,"{%ld}",yip[i]);  /* user id # */
171      FPRINTF(stderr, "Integration will be terminated at the end of ");      name = WriteInstanceNameString(in,g_solvinst_cur);
172      FPRINTF(stderr, "the current step.\n");      FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));
173      return 0;      ascfree(name);
174    }    }
175    if (status.time_limit_exceeded) {    FPRINTF(fp,BCOLSFMT,"indvar");
176      FPRINTF(stderr, "The time limit was ");    for (i=0; i < len; i++) {
177      FPRINTF(stderr, "exceeded while calculating derivatives.\n");      FPRINTF(fp,BCOLIFMT,yip[i]);
178      FPRINTF(stderr, "Integration will be terminated at ");    }
179      FPRINTF(stderr, "the end of the current step.\n");    FPRINTF(fp,"\n");
180      return 0;    for (i=0; i <= len; i++) {
181    }      FPRINTF(fp,BCOLHFMT);
182    if (status.iteration_limit_exceeded) {    }
183      FPRINTF(stderr, "The iteration limit was ");    FPRINTF(fp,"\n");
184      FPRINTF(stderr, "exceeded while calculating derivatives.\n");  }
185      FPRINTF(stderr, "Integration will be terminated at ");  /********************************************************************/
186      FPRINTF(stderr, "the end of the current step.\n");  void Asc_IntegPrintObsHeader(FILE *fp, IntegratorSystem *blsys)
187      return 0;  {
188    }    long i,len, *obsip;
189    if (status.panic) {    char *name;
190      FPRINTF(stderr, "The user patience limit was ");    struct Instance *in;
191      FPRINTF(stderr, "exceeded while calculating derivatives.\n");    int si;
192      FPRINTF(stderr, "Integration will be terminated at ");    if (fp==NULL) {
193      FPRINTF(stderr, "the end of the current step.\n");      return;
194      return 0;    }
195    }    if (blsys==NULL) {
196    return 0;      ERROR_REPORTER_HERE(ASC_PROG_ERR,"called without data");
197  }      return;
198      }
199  /********************************************************************/    if (blsys->n_obs == 0) {
200  static FILE *l_y_file = NULL;      return;
201  static FILE *l_obs_file = NULL;    }
202  static char *l_y_filename = NULL;    if (blsys->obs == NULL) {
203  static char *l_obs_filename = NULL;      ERROR_REPORTER_HERE(ASC_PROG_ERR,"called with NULL data");
204  static int l_print_option = 1; /* default si */      return;
205  static int l_print_fmt = 0; /* default variable */    }
206      len = blsys->n_obs;
207  FILE *Asc_IntegOpenYFile(void)    obsip = blsys->obs_id;
208  {    si = l_print_option;
209    if (l_y_filename==NULL) {    FPRINTF(fp,"Observations: (user index) (name) (units)\n");
210      return NULL;    /* output indep var name */
211    }    /* output obs var names */
212    l_y_file = fopen(l_y_filename,"a+");    in = var_instance(blsys->x);
213      FPRINTF(fp,"{indvar}");  /* indep id name */
214    if (l_y_file==NULL) {    name = WriteInstanceNameString(in,g_solvinst_cur);
215      FPRINTF(ASCERR,    FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));
216        "WARNING: (integrate) Unable to open\n\t%s\nfor state output log.\n",    ascfree(name);
217        l_y_filename);    for (i=0; i< len; i++) {
218    } else {      in = var_instance(blsys->obs[i]);
219      time_t t;      FPRINTF(fp,"{%ld}",obsip[i]);  /* user id # */
220      t = time((time_t *)NULL);      name = WriteInstanceNameString(in,g_solvinst_cur);
221      FPRINTF(l_y_file,"DATASET %s", asctime(localtime(&t)));      FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));
222      FFLUSH(l_y_file);      ascfree(name);
223    }    }
224    return l_y_file;    FPRINTF(fp,BCOLSFMT,"indvar");
225  }    for (i=0; i < len; i++) {
226        FPRINTF(fp,BCOLIFMT,obsip[i]);
227  FILE *Asc_IntegOpenObsFile(void)    }
228  {    FPRINTF(fp,"\n");
229    if (l_obs_filename==NULL) {    for (i=0; i <= len; i++) {
230      return NULL;      FPRINTF(fp,BCOLHFMT);
231    }    }
232    l_obs_file = fopen(l_obs_filename,"a+");    FPRINTF(fp,"\n");
233    if (l_obs_file==NULL) {  }
234      FPRINTF(ASCERR,  
235        "WARNING: (integrate) Unable to open\n\t%s\nfor observation log.\n",  /********************************************************************/
236        l_obs_filename);  void Asc_IntegPrintYLine(FILE *fp, IntegratorSystem *blsys)
237    } else {  {
238      time_t t;    long i,len;
239      t = time((time_t *)NULL);    struct var_variable **vp;
240      FPRINTF(l_obs_file,"DATASET %s", asctime(localtime(&t)));    int si;
241      FFLUSH(l_obs_file);    if (fp==NULL) {
242    }      return;
243    return l_obs_file;    }
244  }    if (blsys==NULL) {
245        FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintYLine: called w/o data\n");
246  FILE *Asc_IntegGetYFile(void)      return;
247  {    }
248    return l_y_file;    if (blsys->n_y == 0) {
249  }      return;
250      }
251  FILE *Asc_IntegGetObsFile(void)    if (blsys->y == NULL) {
252  {      FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintYHeader: called w/NULL data\n");
253    return l_obs_file;      return;
254  }    }
255      vp = blsys->y;
256  void Asc_IntegReleaseYFile(void)    len = blsys->n_y;
257  {    si = l_print_option;
258    l_y_file = NULL;    FPRINTF(fp,BCOLNFMT,Asc_UnitlessValue(var_instance(blsys->x),si));
259  }    for (i=0; i < len; i++) {
260  void Asc_IntegReleaseObsFile(void)      FPRINTF(fp,BCOLNFMT, Asc_UnitlessValue(var_instance(vp[i]),si));
261  {    }
262    l_obs_file = NULL;    FPRINTF(fp,"\n");
263  }  }
264  /********************************************************************/  
265  /* string column layout: 26 char columns  if l_print_fmt else variable */  void Asc_IntegPrintObsLine(FILE *fp, IntegratorSystem *blsys)
266  /* numeric format: leading space plus characters from Unit*Value */  {
267  /* index format: left padded long int */    long i,len;
268  /* headline format space plus - s */    struct var_variable **vp;
269  #define BCOLSFMT (l_print_fmt ? "%-26s" : "\t%s")    int si;
270  #define BCOLNFMT (l_print_fmt ? " %-25s" : "\t%s")    if (fp==NULL) {
271  #define BCOLIFMT (l_print_fmt ? " %25ld" : "\t%ld")      return;
272  #define BCOLHFMT (l_print_fmt ? " -------------------------" : "\t---")    }
273      if (blsys==NULL) {
274  void Asc_IntegPrintYHeader(FILE *fp, struct Integ_system_t *blsys)      FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintObsLine: called w/o data\n");
275  {      return;
276    long i,len, *yip;    }
277    char *name;    if (blsys->n_obs == 0) {
278    struct Instance *in;      return;
279    int si;    }
280    if (fp==NULL) {    if (blsys->obs == NULL) {
281      return;      FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintObsHeader: called w/NULL data\n");
282    }      return;
283    if (blsys==NULL) {    }
284      FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintYHeader: called w/o data\n");    vp = blsys->obs;
285      return;    len = blsys->n_obs;
286    }    si = l_print_option;
287    if (blsys->n_y == 0) {    FPRINTF(fp,BCOLNFMT,Asc_UnitlessValue(var_instance(blsys->x),si));
288      return;    for (i=0; i < len; i++) {
289    }      FPRINTF(fp,BCOLNFMT, Asc_UnitlessValue(var_instance(vp[i]),si));
290    if (blsys->y == NULL) {    }
291      FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintYHeader: called w/NULL data\n");    FPRINTF(fp,"\n");
292      return;  }
293    }  
294    len = blsys->n_y;  
295    yip = blsys->y_id;  /*---------------------------------------------*/
296    si = l_print_option;  
297    /* output indep var name */  int Asc_IntegSetYFileCmd(ClientData cdata, Tcl_Interp *interp,
298    /* output dep var names */                        int argc, CONST84 char *argv[])
299    FPRINTF(fp,"States: (user index) (name) (units)\n");  {
300    in = var_instance(blsys->x);    size_t len;
301    FPRINTF(fp,"{indvar}");  /* indep id name */  
302    name = WriteInstanceNameString(in,g_solvinst_cur);    UNUSED_PARAMETER(cdata);
303    FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));  
304    ascfree(name);    if ( argc != 2 ) {
305    for (i=0; i< len; i++) {      FPRINTF(ASCERR, "integrate_set_y_file: called without filename.\n");
306      in = var_instance(blsys->y[i]);      Tcl_SetResult(interp,
307      FPRINTF(fp,"{%ld}",yip[i]);  /* user id # */                    "integrate_set_y_file <filename,""> called without arg.",
308      name = WriteInstanceNameString(in,g_solvinst_cur);                    TCL_STATIC);
309      FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));      return TCL_ERROR;
310      ascfree(name);    }
311    }    if (len >0 ) {
312    FPRINTF(fp,BCOLSFMT,"indvar");      l_y_filename = Asc_MakeInitString((int)len);
313    for (i=0; i < len; i++) {      sprintf(l_y_filename,"%s",argv[1]);
314      FPRINTF(fp,BCOLIFMT,yip[i]);    } else {
315    }      l_y_filename = NULL;
316    FPRINTF(fp,"\n");    }
317    for (i=0; i <= len; i++) {    return TCL_OK;
318      FPRINTF(fp,BCOLHFMT);  }
319    }  
320    FPRINTF(fp,"\n");  /*---------------------------------------------*/
321  }  
322  /********************************************************************/  int Asc_IntegSetObsFileCmd(ClientData cdata, Tcl_Interp *interp,
323  void Asc_IntegPrintObsHeader(FILE *fp, struct Integ_system_t *blsys)                          int argc, CONST84 char *argv[])
324  {  {
325    long i,len, *obsip;    size_t len;
326    char *name;  
327    struct Instance *in;    UNUSED_PARAMETER(cdata);
328    int si;  
329    if (fp==NULL) {    if ( argc != 2 ) {
330      return;      FPRINTF(ASCERR, "integrate_set_obs_file: called without filename.\n");
331    }      Tcl_SetResult(interp,
332    if (blsys==NULL) {                    "integrate_set_obs_file <filename,""> called without arg.",
333      FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintObsHeader: called w/o data\n");                    TCL_STATIC);
334      return;      return TCL_ERROR;
335    }    }
336    if (blsys->n_obs == 0) {    if (l_obs_filename != NULL) {
337      return;      ascfree(l_obs_filename);
338    }    }
339    if (blsys->obs == NULL) {    len = strlen(argv[1]);
340      FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintObsHeader: called w/NULL data\n");    if (len >0 ) {
341      return;      l_obs_filename = Asc_MakeInitString((int)len);
342    }      sprintf(l_obs_filename,"%s",argv[1]);
343    len = blsys->n_obs;    } else {
344    obsip = blsys->obs_id;      l_obs_filename = NULL;
345    si = l_print_option;    }
346    FPRINTF(fp,"Observations: (user index) (name) (units)\n");    return TCL_OK;
347    /* output indep var name */  }
348    /* output obs var names */  
349    in = var_instance(blsys->x);  /*---------------------------------------------*/
350    FPRINTF(fp,"{indvar}");  /* indep id name */  
351    name = WriteInstanceNameString(in,g_solvinst_cur);  int Asc_IntegSetFileUnitsCmd(ClientData cdata, Tcl_Interp *interp,
352    FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));                              int argc, CONST84 char *argv[])
353    ascfree(name);  {
354    for (i=0; i< len; i++) {    UNUSED_PARAMETER(cdata);
355      in = var_instance(blsys->obs[i]);  
356      FPRINTF(fp,"{%ld}",obsip[i]);  /* user id # */    if ( argc != 2 ) {
357      name = WriteInstanceNameString(in,g_solvinst_cur);      FPRINTF(ASCERR, "integrate_logunits: called without printoption.\n");
358      FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));      Tcl_SetResult(interp,"integrate_logunits <display,si> called without arg.",
359      ascfree(name);                    TCL_STATIC);
360    }      return TCL_ERROR;
361    FPRINTF(fp,BCOLSFMT,"indvar");    }
362    for (i=0; i < len; i++) {    switch (argv[1][0]) {
363      FPRINTF(fp,BCOLIFMT,obsip[i]);    case 's':
364    }      l_print_option = 1;
365    FPRINTF(fp,"\n");      break;
366    for (i=0; i <= len; i++) {    case 'd':
367      FPRINTF(fp,BCOLHFMT);      l_print_option = 0;
368    }      break;
369    FPRINTF(fp,"\n");    default:
370  }      FPRINTF(ASCERR,"integrate_logunits: called with bogus argument.\n");
371        FPRINTF(ASCERR,"logunits remain set to %s.\n",
372  /********************************************************************/        (l_print_option ? "si":"display"));
373  void Asc_IntegPrintYLine(FILE *fp, struct Integ_system_t *blsys)      break;
374  {    }
375    long i,len;    return TCL_OK;
376    struct var_variable **vp;  }
377    int si;  
378    if (fp==NULL) {  /*---------------------------------------------*/
379      return;  
380    }  int Asc_IntegSetFileFormatCmd(ClientData cdata, Tcl_Interp *interp,
381    if (blsys==NULL) {                               int argc, CONST84 char *argv[])
382      FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintYLine: called w/o data\n");  {
383      return;    UNUSED_PARAMETER(cdata);
384    }  
385    if (blsys->n_y == 0) {    if ( argc != 2 ) {
386      return;      FPRINTF(ASCERR, "integrate_logformat called without printoption.\n");
387    }      Tcl_SetResult(interp,
388    if (blsys->y == NULL) {                    "integrate_logformat <fixed,variable> called without arg.",
389      FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintYHeader: called w/NULL data\n");                    TCL_STATIC);
390      return;      return TCL_ERROR;
391    }    }
392    vp = blsys->y;    switch (argv[1][0]) {
393    len = blsys->n_y;    case 'f':
394    si = l_print_option;      l_print_fmt = 1;
395    FPRINTF(fp,BCOLNFMT,Asc_UnitlessValue(var_instance(blsys->x),si));      break;
396    for (i=0; i < len; i++) {    case 'v':
397      FPRINTF(fp,BCOLNFMT, Asc_UnitlessValue(var_instance(vp[i]),si));      l_print_fmt = 0;
398    }      break;
399    FPRINTF(fp,"\n");    default:
400  }      FPRINTF(ASCERR,"integrate_logformat: called with bogus argument.\n");
401        FPRINTF(ASCERR,"logformat remains set to %s.\n",
402  void Asc_IntegPrintObsLine(FILE *fp, struct Integ_system_t *blsys)        (l_print_fmt ? "fixed":"variable"));
403  {      break;
404    long i,len;    }
405    struct var_variable **vp;    return TCL_OK;
406    int si;  }
407    if (fp==NULL) {  
408      return;  
409    }  /*---------------------------------------------------------------
410    if (blsys==NULL) {    DEFINING THE TIMESTEPS FOR INTEGRATION
411      FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintObsLine: called w/o data\n");  */
412      return;  
413    }  int Asc_IntegGetXSamplesCmd(ClientData cdata, Tcl_Interp *interp,
414    if (blsys->n_obs == 0) {                            int argc, CONST84 char *argv[])
415      return;  {
416    }    static char sval[40]; /* buffer long enough to hold a double printed */
417    if (blsys->obs == NULL) {    struct Units *du = NULL;
418      FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintObsHeader: called w/NULL data\n");    const dim_type *dp;
419      return;    long i,len;
420    }    double *uvalues = NULL;
421    vp = blsys->obs;    char *ustring;
422    len = blsys->n_obs;    double *uv;
423    si = l_print_option;    int trydu=0, prec, stat=0;
424    FPRINTF(fp,BCOLNFMT,Asc_UnitlessValue(var_instance(blsys->x),si));    SampleList *samples;
425    for (i=0; i < len; i++) {  
426      FPRINTF(fp,BCOLNFMT, Asc_UnitlessValue(var_instance(vp[i]),si));    UNUSED_PARAMETER(cdata);
427    }  
428    FPRINTF(fp,"\n");    if (( argc < 1 ) || ( argc > 2 )) {
429  }      Tcl_SetResult(interp,
430                      "integrate_get_samples: expected 0 or 1 args [display]",
431  int Asc_IntegSetYFileCmd(ClientData cdata, Tcl_Interp *interp,                    TCL_STATIC);
432                        int argc, CONST84 char *argv[])      return TCL_ERROR;
433  {    }
434    size_t len;    if (argc==2) {
435        trydu = 1;
436    (void)cdata;    /* stop gcc whine about unused parameter */      if( argv[1][0] != 'd') {
437          Tcl_SetResult(interp, "integrate_get_samples: expected display but got ",
438    if ( argc != 2 ) {                      TCL_STATIC);
439      FPRINTF(ASCERR, "integrate_set_y_file: called without filename.\n");        Tcl_AppendResult(interp,argv[1],".",SNULL);
440      Tcl_SetResult(interp,        return TCL_ERROR;
441                    "integrate_set_y_file <filename,""> called without arg.",      }
442                    TCL_STATIC);    }
443      return TCL_ERROR;  
444    }    len = samplelist_length(&l_samplelist);
445    if (l_y_filename != NULL) {    dp = samplelist_dim(&l_samplelist);
446      ascfree(l_y_filename);  
447    }    if (len <1) {
448    len = strlen(argv[1]);      Tcl_SetResult(interp, "{} {}", TCL_STATIC);
449    if (len >0 ) {      return TCL_OK;
450      l_y_filename = Asc_MakeInitString((int)len);    }
451      sprintf(l_y_filename,"%s",argv[1]);  
452    } else {    if (trydu) {
453      l_y_filename = NULL;  
454    }      /* Allocate the space for the retrieved values... */
455    return TCL_OK;      uvalues = ASC_NEW_ARRAY(double,len);
456  }      if (uvalues == NULL) {
457          Tcl_SetResult(interp, "integrate_get_samples: Insufficient memory.",
458  int Asc_IntegSetObsFileCmd(ClientData cdata, Tcl_Interp *interp,                      TCL_STATIC);
459                          int argc, CONST84 char *argv[])        return TCL_ERROR;
460  {      }
461    size_t len;  
462        /* Get the display units at string... */
463    (void)cdata;    /* stop gcc whine about unused parameter */      ustring = Asc_UnitDimString(dp,0);
464    
465    if ( argc != 2 ) {      /* Get the conversion factor... */
466      FPRINTF(ASCERR, "integrate_set_obs_file: called without filename.\n");      du = (struct Units *)LookupUnits(ustring);
467      Tcl_SetResult(interp,      if (du == NULL) {
468                    "integrate_set_obs_file <filename,""> called without arg.",        ERROR_REPORTER_HERE(ASC_PROG_ERR,"LookupUnits failed :-/");
469                    TCL_STATIC);        stat = 1;
470      return TCL_ERROR;      }else{
471    }        stat = 0;
472    if (l_obs_filename != NULL) {        /* fill 'uvalues' with scaled values (in output units) */
473      ascfree(l_obs_filename);        uv = uvalues;
474    }        for (i=0; i < len; i++) {
475    len = strlen(argv[1]);          /* convert to output units */
476    if (len >0 ) {          stat = Asc_UnitConvert(du,samplelist_get(&l_samplelist,i),uv,1);
477      l_obs_filename = Asc_MakeInitString((int)len);          if (stat) {
478      sprintf(l_obs_filename,"%s",argv[1]);            /* any problems, just stop */
479    } else {            break;
480      l_obs_filename = NULL;          }
481    }          uv++;
482    return TCL_OK;        }
483  }      }
484        if (stat) {
485  int Asc_IntegSetFileUnitsCmd(ClientData cdata, Tcl_Interp *interp,        /* there was a problem, so free the allocated space */
486                              int argc, CONST84 char *argv[])        ascfree(uvalues);
487  {      }
488    (void)cdata;    /* stop gcc whine about unused parameter */    }
489    
490    if ( argc != 2 ) {    /* give Tcl the units string */
491      FPRINTF(ASCERR, "integrate_logunits: called without printoption.\n");    Tcl_AppendResult(interp,"{",ustring,"} {",SNULL);
492      Tcl_SetResult(interp,"integrate_logunits <display,si> called without arg.",  
493                    TCL_STATIC);    /* work out what precision we want */
494      return TCL_ERROR;    prec = Asc_UnitGetCPrec();
495    }  
496    switch (argv[1][0]) {    len--; /* last one is a special case...? */
497    case 's':    if (!trydu || stat){
498      l_print_option = 1;      /* no unit conversion, or failed unit conversion: use the raw values */
499      break;      for(i=0;i<len;i++){
500    case 'd':          sprintf(sval,"%.*g ",prec,samplelist_get(&l_samplelist,i));
501      l_print_option = 0;          Tcl_AppendResult(interp,sval,SNULL);
502      break;      }
503    default:    }else{
504      FPRINTF(ASCERR,"integrate_logunits: called with bogus argument.\n");      for(i=0;i<len;i++){
505      FPRINTF(ASCERR,"logunits remain set to %s.\n",          sprintf(sval,"%.*g ",prec,uvalues[i]);
506        (l_print_option ? "si":"display"));          Tcl_AppendResult(interp,sval,SNULL);
507      break;      }
508    }      ascfree(uvalues);
509    return TCL_OK;    }
510  }    for (i=0; i<len; i++) {
511        /* print number and a blank */
512  int Asc_IntegSetFileFormatCmd(ClientData cdata, Tcl_Interp *interp,      sprintf(sval,"%.*g ",prec,uvalues[i]);
513                               int argc, CONST84 char *argv[])      Tcl_AppendResult(interp,sval,SNULL);
514  {    }
515    (void)cdata;    /* stop gcc whine about unused parameter */  
516      sprintf(sval,"%.*g",prec,uvalues[len]);
517    if ( argc != 2 ) {    Tcl_AppendResult(interp,sval,"}",SNULL);
518      FPRINTF(ASCERR, "integrate_logformat called without printoption.\n");  
519      Tcl_SetResult(interp,    if (trydu && !stat) {
520                    "integrate_logformat <fixed,variable> called without arg.",      ascfree(uvalues);
521                    TCL_STATIC);    }
522      return TCL_ERROR;    return TCL_OK;
523    }  }
524    switch (argv[1][0]) {  
525    case 'f':  int Asc_IntegSetXSamplesCmd(ClientData cdata, Tcl_Interp *interp,
526      l_print_fmt = 1;                            int argc, CONST84 char *argv[])
527      break;  {
528    case 'v':    struct Units *du = NULL;
529      l_print_fmt = 0;    const dim_type *dp;
530      break;    dim_type *mydp;
531    default:    long i,len;
532      FPRINTF(ASCERR,"integrate_logformat: called with bogus argument.\n");    double *uvalues = NULL;
533      FPRINTF(ASCERR,"logformat remains set to %s.\n",    double *uv;
534        (l_print_fmt ? "fixed":"variable"));    int stat;
535      break;  
536    }    UNUSED_PARAMETER(cdata);
537    return TCL_OK;  
538  }    if (argc == 1) {
539        samplelist_assign(&l_samplelist,0L,NULL,NULL);
540        return TCL_OK;
541  /********************************************************************/    }
542  /* stuff for handling a user defined list of step outside the ascend  
543     model */    if (argc <4) {
544        Tcl_SetResult(interp,
545  int Asc_IntegSetXSamples(long n, double *values, dim_type *d)                    "Syntax: integrate_set_samples"
546  {                    " <units> <value [value...] value>",
547    /* trash the old stuff, regardless.*/                    TCL_STATIC);
548    g_intg_samples.ns = 0L;      FPRINTF(ASCERR,"ERROR: integrate_set_samples needs at least 3 args.");
549    if (g_intg_samples.x != NULL) {      return TCL_ERROR;
550      ascfree(g_intg_samples.x);    }
551      g_intg_samples.x = NULL;  
552    }    du = (struct Units *)LookupUnits(argv[1]);
553    if (g_intg_samples.d != NULL) {    if (du == NULL) {
554      ascfree(g_intg_samples.d);      Tcl_SetResult(interp, "integrate_set_samples: first arg not valid units.",
555      g_intg_samples.d = NULL;                    TCL_STATIC);
556    }      return TCL_ERROR;
557    /* if was a reset, return now */    }
558    if (n <1 || values==NULL) {    dp = (const dim_type *)UnitsDimensions(du);
559      return 0;  #if INTDEBUG
560    }    FPRINTF(ASCERR,"user dimen looks like\n");
561    /* store d given or copy and store WildDimension */    PrintDimen(ASCERR,dp);
562    if (d != NULL) {  #endif
563      g_intg_samples.d = d;    mydp = ASC_NEW(dim_type); /* we are giving away */
564  #if INTDEBUG    if (mydp == NULL) {
565      FPRINTF(ASCERR,"sample received dimen\n");      Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory",
566      PrintDimen(ASCERR,g_intg_samples.d);                    TCL_STATIC);
567  #endif      return TCL_ERROR;
568    } else {    }
569      g_intg_samples.d = (dim_type *)ascmalloc(sizeof(dim_type));    CopyDimensions(dp,mydp);
570      if (g_intg_samples.d == NULL) {  #if INTDEBUG
571        FPRINTF(ASCERR,"ERROR: (Asc_IntegSetXSamples) Insufficient memory.\n");    FPRINTF(ASCERR,"copy of user dimen looks like\n");
572        return 1;    PrintDimen(ASCERR,mydp);
573      }  #endif
574      CopyDimensions(WildDimension(),g_intg_samples.d);  
575  #if INTDEBUG    len = argc -2;
576      FPRINTF(ASCERR,"copy of wild dimen looks like\n");    uvalues = ASC_NEW_ARRAY(double,len); /* we are giving away */
577      PrintDimen(ASCERR,g_intg_samples.d);    if (uvalues==NULL) {
578  #endif      Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory",
579    }                    TCL_STATIC);
580    g_intg_samples.ns = n;      ascfree(mydp);
581    g_intg_samples.x = values;      return TCL_ERROR;
582    return 0;    }
583  }    stat = 0;
584      uv = uvalues;
585  int Asc_IntegGetMaxSteps(struct Integ_system_t *sys)    for (i=0; i<len; i++) {
586  {      if(Tcl_GetDouble(interp,argv[i+2],uv)!=TCL_OK) {
587    return sys->maxsteps;        stat = 1;
588  }        break;
589        }
590  double Asc_IntegGetStepMax(struct Integ_system_t *sys)      stat = Asc_UnitConvert(du,*uv,uv,0);
591  {      if (stat) {
592    return sys->stepmax;        break;
593  }      }
594        uv++;
595  double Asc_IntegGetStepMin(struct Integ_system_t *sys)    }
596  {    Tcl_ResetResult(interp);
597    return sys->stepmin;    if (stat) {
598  }      Tcl_SetResult(interp, "integrate_set_samples: Invalid value given. (",
599                      TCL_STATIC);
600  double Asc_IntegGetStepZero(struct Integ_system_t *sys)      Tcl_AppendResult(interp,argv[i+2],")",SNULL);
601  {      ascfree(uvalues);
602    return sys->stepzero;      ascfree(mydp);
603  }      return TCL_ERROR;
604      }
605  long Asc_IntegGetNSamples()    if(!samplelist_assign(&l_samplelist,len,uvalues,mydp)){
606  {      Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory.",
607    return g_intg_samples.ns;                    TCL_STATIC);
608  }      ascfree(uvalues);
609        ascfree(mydp);
610  double Asc_IntegGetXSamplei(long i)      return TCL_ERROR;
611  {    }
612    if (i > -1 && i < g_intg_samples.ns) {    return TCL_OK;
613      return g_intg_samples.x[i];  }
614    } else {  
615      FPRINTF(ASCERR,  /*----------------------------------------------------------------
616              "WARNING: (Asc_IntegGetXSamplei) Undefined xsample %ld."    FUNCTIONS THAT QUERY THE INTEGRATOR/SOLVER
617              " Returning 0.0.\n",  */
618              i);  
619      return (double)0.0;  /**
620    }      Tcl/Tk interface function: is the specified instance integrable.
621  }  
622        There is a problem with this part of the interface. The
623  void Asc_IntegSetXSamplei(long i,double xi)      'integrator_isintegrable' function is being called before the
624  {      'integrator_analyse' step, which means that it really doesn't have all the
625    if (i > -1 && i < g_intg_samples.ns) {      information that it needs to be able to say. I've reworked it so that
626      g_intg_samples.x[i] = xi;      you would expect to get errors back from the 'integrator_solve' step if
627    } else {      there is any problem.
628      FPRINTF(ASCERR,  
629        "WARNING: (Asc_IntegSetXSamplei) Undefined xsample %ld. Ignored.\n",      All this function now does is to check that the instance is OK, and to
630        i);      check that a valid integrator engine is specified.
631    }  */
632  }  int Asc_IntegInstIntegrableCmd(ClientData cdata,Tcl_Interp *interp,
633                                 int argc, CONST84 char *argv[])
634  dim_type *Asc_IntegGetXSamplesDim(void)  {
635  {    struct Instance *i=NULL;
636    return g_intg_samples.d;    IntegratorEngine integrator=INTEG_UNKNOWN;
637  }    int result=0;         /* 0 = FALSE; 1 = TRUE */
638  /********************************************************************/  
639  /* callbacks for manipulating a user defined list of steps */    UNUSED_PARAMETER(cdata);
640  int Asc_IntegGetXSamplesCmd(ClientData cdata, Tcl_Interp *interp,  
641                            int argc, CONST84 char *argv[])    if ( argc != 3 ) {
642  {      Tcl_SetResult(interp, "integrate_able <solver,current,search> <lsode>",
643    static char sval[40]; /* buffer long enough to hold a double printed */                    TCL_STATIC);
644    struct Units *du = NULL;      return TCL_ERROR;
645    dim_type *dp;    }
646    long i,len;  
647    double *uvalues = NULL;    if (strncmp(argv[1],"solver",3)==0) {
648    char *ustring;      i = g_solvinst_cur;
649    double *uv;    } else {
650    int trydu=0, prec, stat=0;      if (strncmp(argv[1],"search",3)==0) {
651          i = g_search_inst;
652    (void)cdata;    /* stop gcc whine about unused parameter */      } else {
653          if (strncmp(argv[1],"current",3)==0) {
654    if (( argc < 1 ) || ( argc > 2 )) {          i = g_curinst;
655      Tcl_SetResult(interp,        } else {
656                    "integrate_get_samples: expected 0 or 1 args [display]",          Tcl_SetResult(interp,
657                    TCL_STATIC);                        "integrate_able: arg 1 is current, search, or solver",
658      return TCL_ERROR;                        TCL_STATIC);
659    }          return TCL_ERROR;
660    if (argc==2) {        }
661      trydu = 1;      }
662      if( argv[1][0] != 'd') {    }
663        Tcl_SetResult(interp, "integrate_get_samples: expected display but got ",  
664                      TCL_STATIC);    if (!i) {
665        Tcl_AppendResult(interp,argv[1],".",SNULL);      Tcl_SetResult(interp, "0", TCL_STATIC);
666        return TCL_ERROR;      FPRINTF(ASCERR,"NULL instance sent to integrate_able.\n");
667      }      return TCL_OK;
668    }    }
669    len = g_intg_samples.ns;  
670    dp = g_intg_samples.d;    integrator = INTEG_UNKNOWN;
671    if (len <1) {  
672      Tcl_SetResult(interp, "{} {}", TCL_STATIC);    if (strncmp(argv[2],"blsode",3)==0) {
673      return TCL_OK;      integrator = INTEG_LSODE;
674    }    }else if (strncmp(argv[2],"ida",3)==0) {
675        integrator = INTEG_IDA;
676    if (trydu) {    }
677      uvalues = (double *)ascmalloc(sizeof(double)*len);  
678      if (uvalues == NULL) {    result = (integrator != INTEG_UNKNOWN);
679        Tcl_SetResult(interp, "integrate_get_samples: Insufficient memory.",    if (result) {
680                      TCL_STATIC);      Tcl_SetResult(interp, "1", TCL_STATIC);
681        return TCL_ERROR;    } else {
682      }      Tcl_SetResult(interp, "0", TCL_STATIC);
683      ustring = Asc_UnitDimString(dp,0); /* get display unit string */    }
684      du = (struct Units *)LookupUnits(ustring);    return TCL_OK;
685      if (du == NULL) {  }
686        /*  a very bizarre thing to happen,  
687         *  since Asc_UnitDimString just made it  /*-----------------------------------*/
688         */  
689        stat = 1;  /**
690      } else {      Set up the Integrator.
691        stat = 0;  
692        uv = uvalues;      switches (in Tcl/Tk)
693        for (i=0; i < len; i++) {          -engine $name
694          stat = Asc_UnitConvert(du,g_intg_samples.x[i],uv,1);          -i0 $stepindex
695          if (stat) {          -i1 $stepindex
696            break;          -dt0 $initstepsize
697          }          -dtmin $minstep
698          uv++;          -dtmax $maxstep
699        }          -moststeps $moststeps
700      }  */
701      if (stat) {  int Asc_IntegSetupCmd(ClientData cdata,Tcl_Interp *interp,
702        ascfree(uvalues);                     int argc, CONST84 char *argv[])
703      }  {
704    }    IntegratorEngine integrator = INTEG_UNKNOWN;
705    /* if display not wanted or display -> conversion error */    char buf[MAXIMUM_NUMERIC_LENGTH];         /* string to hold integer */
706    if (!trydu || stat) {    CONST84 char *engine = NULL;
707      uvalues = g_intg_samples.x;    int result = 0;         /* 0 = FALSE; 1 = TRUE */
708      ustring = Asc_UnitDimString(dp,1); /* get si unit string */    long i0=(-1), i1=(-1);
709    }    int ifound = 0;
710      int k;
711    Tcl_AppendResult(interp,"{",ustring,"} {",SNULL);    int moststeps=0;
712    prec = Asc_UnitGetCPrec();    double dt0=0, dtmin=0, dtmax=0;
713    len--; /* last one is a special case */    CONST84 char *cdt0=NULL, *cdtmin=NULL, *cdtmax=NULL, *cmoststeps=NULL,
714    for (i=0; i<len; i++) {         *ci0=NULL, *ci1=NULL;
715      /* print number and a blank */    IntegratorReporter *reporter;
716      sprintf(sval,"%.*g ",prec,uvalues[i]);    IntegratorSystem *blsys;
717      Tcl_AppendResult(interp,sval,SNULL);  
718    }    UNUSED_PARAMETER(cdata);
719    sprintf(sval,"%.*g",prec,uvalues[len]);  
720    Tcl_AppendResult(interp,sval,"}",SNULL);    k = 1;
721    if (trydu && !stat) {    while (k < (argc-1)) { /* arguments come in pairs */
722      ascfree(uvalues);      if (strcmp(argv[k],"-engine")==0) {
723    }        engine = argv[k+1];
724    return TCL_OK;        k+=2;
725  }        continue;
726        }
727  int Asc_IntegSetXSamplesCmd(ClientData cdata, Tcl_Interp *interp,      if (strcmp(argv[k],"-i1")==0) {
728                            int argc, CONST84 char *argv[])        ci1 = argv[k+1];
729  {        k+=2;
730    struct Units *du = NULL;        continue;
731    dim_type *dp;      }
732    dim_type *mydp;      if (strcmp(argv[k],"-i0")==0) {
733    long i,len;        ci0 = argv[k+1];
734    double *uvalues = NULL;        k+=2;
735    double *uv;        continue;
736    int stat;      }
737        if (strcmp(argv[k],"-moststeps")==0) {
738    (void)cdata;    /* stop gcc whine about unused parameter */        cmoststeps = argv[k+1];
739          k+=2;
740    if (argc == 1) {        continue;
741      Asc_IntegSetXSamples(0L,NULL,NULL);      }
742      return TCL_OK;      if (strcmp(argv[k],"-dtmax")==0) {
743    }        cdtmax = argv[k+1];
744          k+=2;
745    if (argc <4) {        continue;
746      Tcl_SetResult(interp,      }
747                    "Syntax: integrate_set_samples"      if (strcmp(argv[k],"-dtmin")==0) {
748                    " <units> <value [value...] value>",        cdtmin = argv[k+1];
749                    TCL_STATIC);        k+=2;
750      FPRINTF(ASCERR,"ERROR: integrate_set_samples needs at least 3 args.");        continue;
751      return TCL_ERROR;      }
752    }      if (strcmp(argv[k],"-dt0")==0) {
753          cdt0 = argv[k+1];
754    du = (struct Units *)LookupUnits(argv[1]);        k+=2;
755    if (du == NULL) {        continue;
756      Tcl_SetResult(interp, "integrate_set_samples: first arg not valid units.",      }
757                    TCL_STATIC);      Tcl_AppendResult(interp,argv[0],": unrecognized option ",
758      return TCL_ERROR;                       argv[k],".",SNULL);
759    }      return TCL_ERROR;
760    dp = (dim_type *)UnitsDimensions(du);    }
761  #if INTDEBUG  
762    FPRINTF(ASCERR,"user dimen looks like\n");    if (engine != NULL && strncmp(engine,"BLSODE",3)==0) {
763    PrintDimen(ASCERR,dp);      integrator = INTEG_LSODE;
764  #endif      ifound=1;
765    mydp = (dim_type *)ascmalloc(sizeof(dim_type)); /* we are giving away */    }
766    if (mydp == NULL) {  
767      Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory",    if (engine != NULL && strncmp(engine,"IDA",3)==0) {
768                    TCL_STATIC);      integrator = INTEG_IDA;
769      return TCL_ERROR;      ifound=1;
770    }    }
771    CopyDimensions(dp,mydp);  
772  #if INTDEBUG    if (!ifound) {
773    FPRINTF(ASCERR,"copy of user dimen looks like\n");      Tcl_SetResult(interp, "Unsupported integrator", TCL_STATIC);
774    PrintDimen(ASCERR,mydp);      Tcl_AppendResult(interp," ",engine,SNULL);
775  #endif      return TCL_ERROR;
776      }
777    len = argc -2;    if (ci0 != NULL && ci1 != NULL) {
778    uvalues = (double *)ascmalloc(sizeof(double) * len); /* we are giving away */      /* get i0, i1 if both supplied. */
779    if (uvalues==NULL) {      long i;
780      Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory",      if (Tcl_ExprLong(interp,ci0,&i)==TCL_ERROR|| i<0) {
781                    TCL_STATIC);        Tcl_ResetResult(interp);
782      ascfree(mydp);        Tcl_SetResult(interp, "integrate_setup: index i0 invalid", TCL_STATIC);
783      return TCL_ERROR;        return TCL_ERROR;
784    }      }
785    stat = 0;      i0=i;
786    uv = uvalues;      if (Tcl_ExprLong(interp,ci1,&i)==TCL_ERROR|| i<i0) {
787    for (i=0; i<len; i++) {        Tcl_ResetResult(interp);
788      if(Tcl_GetDouble(interp,argv[i+2],uv)!=TCL_OK) {        Tcl_SetResult(interp, "integrate_setup: index i1 invalid", TCL_STATIC);
789        stat = 1;        return TCL_ERROR;
790        break;      }
791      }      i1=i;
792      stat = Asc_UnitConvert(du,*uv,uv,0);    }
793      if (stat) {    if (cdt0 != NULL) {
794        break;      if (Tcl_GetDouble(interp,cdt0,&dt0) != TCL_OK) {
795      }        Tcl_ResetResult(interp);
796      uv++;        Tcl_AppendResult(interp, "integrate_setup: initial step length invalid",
797    }                         " (",cdt0,")", SNULL);
798    Tcl_ResetResult(interp);        return TCL_ERROR;
799    if (stat) {      }
800      Tcl_SetResult(interp, "integrate_set_samples: Invalid value given. (",    }
801                    TCL_STATIC);    if (cdtmin != NULL) {
802      Tcl_AppendResult(interp,argv[i+2],")",SNULL);      if (Tcl_GetDouble(interp,cdtmin,&dtmin) != TCL_OK || dtmin < 0) {
803      ascfree(uvalues);        Tcl_ResetResult(interp);
804      ascfree(mydp);        Tcl_AppendResult(interp, "integrate_setup: minimum step length invalid",
805      return TCL_ERROR;                         " (",cdtmin,")", SNULL);
806    }        return TCL_ERROR;
807    if (Asc_IntegSetXSamples(len,uvalues,mydp)) {      }
808      Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory.",    }
809                    TCL_STATIC);    if (cdtmax != NULL) {
810      ascfree(uvalues);      if (Tcl_GetDouble(interp,cdtmax,&dtmax) != TCL_OK || dtmax < dtmin) {
811      ascfree(mydp);        Tcl_ResetResult(interp);
812      return TCL_ERROR;        Tcl_AppendResult(interp, "integrate_setup: maximum step length invalid",
813    }                         " (",cdtmax,")", SNULL);
814    return TCL_OK;        return TCL_ERROR;
815  }      }
816      }
817  /********************************************************************/    if (cmoststeps != NULL) {
818        if (Tcl_GetInt(interp,cmoststeps,&moststeps) != TCL_OK || moststeps < 0) {
819  int Asc_IntegInstIntegrable(struct Instance *inst,        Tcl_ResetResult(interp);
820                            enum Integrator_type integrator)        Tcl_AppendResult(interp, "integrate_setup: maximum internal steps bad",
821  {                         " (",cmoststeps,")", SNULL);
822    switch (integrator) {        return TCL_ERROR;
823    case BLSODE :      }
824      if (inst == NULL) {    }
825        return 0;  
826      }    reporter = Asc_GetIntegReporter();
827      return 1;  
828    case UNKNOWN:    blsys = integrator_new(g_solvsys_cur,g_solvinst_cur);
829      FPRINTF(stderr, "UNKNOWN integrator is not supported.\n");    integrator_set_engine(blsys, integrator);
830      return 0;    integrator_set_reporter(blsys, reporter);
831    default :    integrator_set_samples(blsys,&l_samplelist);
832      FPRINTF(stderr, "ERRONEOUS integrator is not supported.\n");    integrator_set_stepzero(blsys,dt0);
833      return 0;    integrator_set_minstep(blsys,dtmin);
834    }    integrator_set_maxstep(blsys,dtmax);
835  }    integrator_set_maxsubsteps(blsys,moststeps);
836    
837  int Asc_IntegInstIntegrableCmd(ClientData cdata,Tcl_Interp *interp,    result = integrator_analyse(blsys);
838                               int argc, CONST84 char *argv[])  
839  {    /* go and solve it */
840    struct Instance *i=NULL;    integrator_solve(blsys, i0, i1);
841    enum Integrator_type integrator=UNKNOWN;  
842    int result=0;         /* 0 = FALSE; 1 = TRUE */    /* once solution is finished, free whatever we allocated */
843      integrator_free(blsys);
844    (void)cdata;    /* stop gcc whine about unused parameter */  
845      sprintf(buf, "%d", result);
846    if ( argc != 3 ) {    Tcl_SetResult(interp, buf, TCL_VOLATILE);
847      Tcl_SetResult(interp, "integrate_able <solver,current,search> <lsode>",    return TCL_OK;
848                    TCL_STATIC);  }
849      return TCL_ERROR;  
850    }  /********************************************************************/
851    
852    if (strncmp(argv[1],"solver",3)==0) {  int Asc_IntegCleanupCmd(ClientData cdata,Tcl_Interp *interp,
853      i = g_solvinst_cur;                       int argc, CONST84 char *argv[])
854    } else {  {
855      if (strncmp(argv[1],"search",3)==0) {    UNUSED_PARAMETER(cdata);
856        i = g_search_inst;    (void)argv;     /* stop gcc whine about unused parameter */
857      } else {  
858        if (strncmp(argv[1],"current",3)==0) {    if (argc!=1) {
859          i = g_curinst;      Tcl_SetResult(interp, "integrate_cleanup takes no arguments", TCL_STATIC);
860        } else {      return TCL_ERROR;
861          Tcl_SetResult(interp,    }
862                        "integrate_able: arg 1 is current, search, or solver",  
863                        TCL_STATIC);    /* integrator_cleanup(); */
864          return TCL_ERROR;    return TCL_OK;
865        }  }
866      }  
867    }  /*-------------------------------------------------------------------
868      REPORTER FUNCTIONS
869    if (!i) {  */
870      Tcl_SetResult(interp, "0", TCL_STATIC);  
871      FPRINTF(ASCERR,"NULL instance sent to integrate_able.\n");  FILE *integ_y_out;
872      return TCL_OK;  FILE *integ_obs_out;
873    }  
874    IntegratorReporter *Asc_GetIntegReporter(){
875    if (strncmp(argv[2],"blsode",3)==0) {      IntegratorReporter *r;
876      integrator = BLSODE;      r = (IntegratorReporter *)ascmalloc(sizeof(IntegratorReporter));
877    }      r->init = &Asc_IntegReporterInit;
878    /* someday will have an ivp homegrown. in anycase, ivp dis/en-ables btn */      r->write = &Asc_IntegReporterWrite;
879    if (strncmp(argv[2],"ivp",3)==0) {      r->write_obs = &Asc_IntegReporterWriteObs;
880      integrator = UNKNOWN;      r->close = &Asc_IntegReporterClose;
881    }      CONSOLE_DEBUG("CREATED INTEGRATORREPORTER FOR TCL/TK INTERFACE");
882        return r;
883    result = Asc_IntegInstIntegrable(i, integrator);  }
884    if (result) {  
885      Tcl_SetResult(interp, "1", TCL_STATIC);  int Asc_IntegReporterInit(IntegratorSystem *blsys){
886    } else {      CONSOLE_DEBUG("INITIALISING REPORTER");
887      Tcl_SetResult(interp, "0", TCL_STATIC);  
888    }      /* set up output files */
889    return TCL_OK;      integ_y_out = Asc_IntegOpenYFile();
890  }      integ_obs_out = Asc_IntegOpenObsFile();
891    
892  /*** derivative parts **** x *****************************************/      CONSOLE_DEBUG("RELEASING FILES");
893    
894  double Asc_IntegGetDX() {      Asc_IntegReleaseYFile();
895    return var_value(g_intginst_d_x);      Asc_IntegReleaseObsFile();
896  }  
897        CONSOLE_DEBUG("WRITING HEADERS");
898  void Asc_IntegSetDX( double value) {  
899    var_set_value(g_intginst_d_x, value);      /* write headers to yout, obsout and initial points */
900    print_debug("set_d_x = %g\n", value);      Asc_IntegPrintYHeader(integ_y_out,blsys);
901  }      Asc_IntegPrintYLine(integ_y_out,blsys);
902        Asc_IntegPrintObsHeader(integ_obs_out,blsys);
903  /*** derivative parts **** y *****************************************/      Asc_IntegPrintObsLine(integ_obs_out,blsys);
904    }
905  double *Asc_IntegGetDY(double *y) {  
906    long i;  int Asc_IntegReporterWrite(IntegratorSystem *blsys){
907        /* write out a line of stuff */
908    if (y==NULL) {      Asc_IntegPrintYLine(integ_y_out,blsys);
909      y = (double *)asccalloc((g_intginst_d_n_eq+1), sizeof(double));  }
910      /* C y[0]  <==> ascend d.y[1]  <==>  f77 y(1) */  
911    }  int Asc_IntegReporterWriteObs(IntegratorSystem *blsys){
912        Asc_IntegPrintObsLine(integ_obs_out,blsys);
913    for (i=1; i<=g_intginst_d_n_eq; i++) {  }
914      y[i-1] = var_value(D_Y_VAR(i));  
915      print_debug("*get_d_y[%ld] = ", i);  int Asc_IntegReporterClose(IntegratorSystem *blsys){
916      print_debug("%g\n", y[i-1]);      /* close the file streams */
917    }      if (integ_y_out!=NULL) {
918    return y;          fclose(integ_y_out);
919  }      }
920    
921  void Asc_IntegSetDY(double *y) {      if (integ_obs_out!=NULL) {
922    long i;          fclose(integ_obs_out);
923        }
924    /* C y[0]  <==> ascend d.y[1]  <==>  f77 y(1) */  }
925    
   for (i=1; i<=g_intginst_d_n_eq; i++) {  
     var_set_value(D_Y_VAR(i),y[i-1]);  
     print_debug("*set_d_y[%ld] = ", i);  
     print_debug("%g\n", y[i-1]);  
   }  
 }  
   
 /*** derivative parts **** dydx *****************************************/  
   
 double *Asc_IntegGetDDydx(double *dydx) {  
   long i;  
   
   if (dydx==NULL) {  
     dydx = (double *)asccalloc((g_intginst_d_n_eq+1), sizeof(double));  
     /* C dydx[0]  <==> ascend d.dydx[1]  <==>  f77 ydot(1) */  
   }  
   
   for (i=1; i<=g_intginst_d_n_eq; i++) {  
     dydx[i-1] = var_value(D_DYDX_VAR(i));  
     print_debug("*get_d_dydx[%ld] = ", i);  
     print_debug("%g\n", dydx[i-1]);  
   }  
   return dydx;  
 }  
   
 void Asc_IntegSetDDydx(double *dydx) {  
   long i;  
   
   /* C dydx[0]  <==> ascend d.dydx[1]  <==>  f77 ydot(1) */  
   
   for (i=1; i<=g_intginst_d_n_eq; i++) {  
     var_set_value(D_DYDX_VAR(i),dydx[i-1]);  
     print_debug("*set_d_dydx[%ld] = ", i);  
     print_debug("%g\n", dydx[i-1]);  
   }  
 }  
   
 /**** derivative parts * d.obs ******************************************/  
   
 /*  
    This function takes the inst in the solver and returns the vector of  
    observation variables that are located in the submodel d.obs array.  
 */  
 double *Asc_IntegGetDObs(double *obsi) {  
   long i;  
   
   if (obsi==NULL) {  
     obsi = (double *)asccalloc((g_intginst_d_n_obs+1),sizeof(double));  
   }  
   
   /* C obsi[0]  <==> ascend d.obs[1] */  
   
   for (i=1; i<=g_intginst_d_n_obs; i++) {  
     obsi[i-1] = var_value(D_OBS_VAR(i));  
     print_debug("*get_d_obs[%ld] = ", i);  
     print_debug("%g\n", obsi[i-1]);  
   }  
   return obsi;  
 }  
   
 /* All the KEEPIP stuff used to be here. its gone to cvs land */  
   
 /*  
  *  takes the type of integrator and start and finish index and calls the  
  *  appropriate integrator  
  */  
 static void Integ_Solve( enum Integrator_type integrator,  
                          int start_index, int finish_index,  
                          struct Integ_system_t *blsys) {  
   switch (integrator) {  
   case BLSODE:  
     Asc_BLsodeIntegrate(g_intgsys_cur,start_index, finish_index, blsys);  
     return;  
   default:  
     FPRINTF(stderr, "The requested integrator is ");  
     FPRINTF(stderr, "not currently available\n.");  
     return;  
   }  
 }  
   
 static double **MakeDenseMatrix(int nrows, int ncols)  
 {  
   int c;  
   double **result;  
   assert(nrows>0);  
   assert(ncols>0);  
   result = (double **)ascmalloc(nrows*sizeof(double *));  
   for (c=0;c<nrows;c++) {  
     result[c] = (double *)asccalloc(ncols,sizeof(double));  
   }  
   return result;  
 }  
   
 static void DestroyDenseMatrix(double **matrix,int nrows)  
 {  
   int c;  
   if (matrix) {  
     for (c=0;c<nrows;c++) {  
       if (matrix[c]) {  
         ascfree((char *)matrix[c]);  
       }  
     }  
     ascfree((char *)matrix);  
   }  
 }  
   
 /*  
  * needs work. Assumes struct Instance* and struct var_variable*  
  * are synonymous, which demonstrates the need for a method to take  
  * an instance and ask the solvers for its global or local index  
  * if var and inst are decoupled.  
  */  
 static  
 int Integ_SetUpDiffs_BLsode(struct Integ_system_t *blsys) {  
   long n_eqns;  
   unsigned long nch,i;  
   
   struct var_variable **vp;  
   int *ip;  
   
   g_intg_diff.n_eqns = n_eqns = blsys->n_y;  
   g_intg_diff.input_indices = (int *)asccalloc(n_eqns, sizeof(int));  
   g_intg_diff.output_indices = (int *)asccalloc(n_eqns, sizeof(int));  
   g_intg_diff.y_vars = NULL;  
   g_intg_diff.ydot_vars = NULL;  
   g_intg_diff.dydx_dx = MakeDenseMatrix(n_eqns,n_eqns);  
   
   
   /*  
    * Let us now process what we consider *inputs* to the problem as  
    * far as ASCEND is concerned; i.e. the state vars or the y_vars's  
    * if you prefer.  
    */  
   nch = n_eqns;  
   vp = g_intg_diff.y_vars =  
     (struct var_variable **)ascmalloc((nch+1)*sizeof(struct var_variable *));  
   ip = g_intg_diff.input_indices;  
   for (i=0;i<nch;i++) {  
     *vp = (struct var_variable *)blsys->y[i];  
     *ip = var_sindex(*vp);  
     vp++;  
     ip++;  
   }  
   *vp = NULL;   /* terminate */  
   
   /*  
    * Let us now go for the outputs, ie the derivative terms.  
    */  
   vp = g_intg_diff.ydot_vars =  
     (struct var_variable **)ascmalloc((nch+1)*sizeof(struct var_variable *));  
   ip = g_intg_diff.output_indices;  
   for (i=0;i<nch;i++) {  
     *vp = (struct var_variable *)blsys->ydot[i];  
     *ip = var_sindex(*vp);  
     vp++;       /* dont assume that a var is synonymous with */  
     ip++;       /* an Instance; that might/will change soon */  
   }  
   *vp = NULL;       /* terminate */  
   
   return 0;  
 }  
   
 /* Build an analytic jacobian for solving the state system */  
 /*  
  * This necessarily ugly piece of code attempts to create a unique  
  * list of relations that explicitly contain the variables in the  
  * given input list. The utility of this information is that we know  
  * exactly which relations must be differentiated, to fill in the  
  * df/dy matrix. If the problem has very few derivative terms, this will  
  * be of great savings. If the problem arose from the discretization of  
  * a pde, then this will be not so useful. The decision wether to use  
  * this function or to simply differentiate the entire relations list  
  * must be done before calling this function. Final Note: the callee  
  * owns the array, but not the array elements.  
  */  
   
 #define AVG_NUM_INCIDENT 4  
   
 /* Returns STATEFLAG child value if we have a correct dynamic variable.  
    If we do not, returns 0. Doesn't check solver_var status.  
    if index pointer is given, it is stuffed too.  
    if return is 0, ignore index.  
    to be correct either type and index are assigned integers or  
    type is -1. */  
 static long DynamicVarInfo(struct var_variable *v,long *index)  
 {  
   struct Instance *c, *d, *i;  
   i = var_instance(v);  
   assert(i!=NULL);  
   c = ChildByChar(i,STATEFLAG);  
   d = ChildByChar(i,STATEINDEX);  
   /* lazy evaluation is important in the following if */  
   if( c == NULL ||  
       d == NULL ||  
       InstanceKind(c) != INTEGER_INST ||  
       InstanceKind(d) != INTEGER_INST ||  
       !AtomAssigned(c) ||  
       (!AtomAssigned(d) && GetIntegerAtomValue(c) != -1L)  
       ) {  
     return 0L;  
   }  
   if (index != NULL) {  
     *index = GetIntegerAtomValue(d);  
   }  
   return GetIntegerAtomValue(c);  
 }  
   
 /* returns the pointer if we have a correct observation variable.  
    If long is passed in, long will have the index value.  
    Vars with UNDEFINED observation flags are 'incorrect.'  
  */  
 static struct var_variable *ObservationVar(struct var_variable *v, long *index)  
 {  
   struct Instance *c,*i;  
   i = var_instance(v);  
   assert(i!=NULL);  
   c = ChildByChar(i,OBSINDEX);  
   if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {  
     return NULL;  
   }  
   if (index != NULL) {  
     *index = GetIntegerAtomValue(c);  
   }  
   return v;  
 }  
   
 /* take the variable instance and set its obsid to index.  
   id must be already assigned at least once. */  
 static void Integ_SetObsId(struct var_variable *v, long index)  
 {  
   struct Instance *c, *i;  
   i = var_instance(v);  
   assert(i!=NULL);  
   c = ChildByChar(i,OBSINDEX);  
   if( c == NULL || InstanceKind(c) != INTEGER_INST || !AtomAssigned(c)) {  
     return;  
   }  
   SetIntegerAtomValue(c,index,0);  
 }  
   
 /*  
  * Sorts the interesting vars into our wonderful little lists.  
  * Dislikes null input intensely.  
  */  
 static void Integ_classify_vars(struct var_variable *var)  
 {  
   struct Integ_var_t *info;  
   long type,index;  
   var_filter_t vfilt;  
   
   assert(var != NULL && var_instance(var)!=NULL );  
   vfilt.matchbits = VAR_SVAR;  
   vfilt.matchvalue = VAR_SVAR;  
   if( var_apply_filter(var,&vfilt) ) {  
     type = DynamicVarInfo(var,&index);  
     if ( type != 0) {  
 #if INTEG_DEBUG  
       var_write_name(g_solvsys_cur,var,ASCERR);  
       FPRINTF(ASCERR," type = %ld\n",type);  
 #endif  
       /* ignore algebraics type == 0 */  
       info = (struct Integ_var_t *)ascmalloc(sizeof(struct Integ_var_t));  
       if (info == NULL) {  
         FPRINTF(ASCERR,  
           "ERROR: (Integ_system_build) Insufficient memory.\n");  
         return;  
       }  
       info->type = type;  
       info->index = index;  
       info->i = var;  
       if (type > 0L) {  
         /* state or derivative */  
         gl_append_ptr(l_isys->dynvars,(POINTER)info);  
         if (type == 1) {  
           l_nstates++;  
         }  
         if (type == 2) {  
           l_nderivs++;  
         }  
       } else {  
         if (type == -1L) {  
           /* independent */  
           gl_append_ptr(l_isys->indepvars,(POINTER)info);  
         } else {  
           /* probably should whine about unknown type here */  
           ascfree(info);  
         }  
       }  
       info = NULL;  
     }  
     if ( ObservationVar(var,&index) != NULL && index > 0L) {  
       info = (struct Integ_var_t *)ascmalloc(sizeof(struct Integ_var_t));  
       if (info == NULL) {  
         FPRINTF(ASCERR,  
           "ERROR: (Integ_system_build) Insufficient memory.\n");  
         return;  
       }  
       info->type = 0L;  
       info->index = index;  
       info->i = var;  
       gl_append_ptr(l_isys->obslist,(POINTER)info);  
       info = NULL;  
     }  
   }  
 }  
   
   
 /* compares observation structs. NULLs should end up at far end. */  
 static int Integ_CmpObs(struct Integ_var_t *v1, struct Integ_var_t *v2)  
 {  
   if (v1 == NULL) {  
     return 1;  
   }  
   if (v2 == NULL) {  
     return -1;  
   }  
   if (v1->index > v2->index) {  
     return 1;  
   }  
   if (v1->index == v2->index) {  
     return 0;  
   }  
   return -1;  
 }  
   
 /*  
   Compares dynamic vars structs. NULLs should end up at far end.  
   List should be sorted primarily by index and then by type, in order  
   of increasing value of both.  
 */  
 static int Integ_CmpDynVars(struct Integ_var_t *v1, struct Integ_var_t *v2)  
 {  
   if (v1 == NULL) {  
     return 1;  
   }  
   if (v2 == NULL) {  
     return -1;  
   }  
   if (v1->index > v2->index) {  
     return 1;  
   }  
   if (v1->index != v2->index) {  
     return -1;  
   }  
   if (v1->type > v2->type) {  
     return 1;  
   }  
   return -1;  
 }  
   
 /* trash nonnull contents of a system_t . does ree the pointer sent.*/  
 static void Integ_system_destroy(struct Integ_system_t *sys)  
 {  
   if (sys==NULL) {  
     return;  
   }  
   if (sys->states != NULL) {  
     gl_destroy(sys->states);  
   }  
   if (sys->derivs != NULL) {  
     gl_destroy(sys->derivs);  
   }  
   if (sys->dynvars != NULL) {  
     gl_free_and_destroy(sys->dynvars);    /* we own the objects in dynvars */  
   }  
   if (sys->obslist != NULL) {  
     gl_free_and_destroy(sys->obslist);    /* and obslist */  
   }  
   if (sys->indepvars != NULL) {  
     gl_free_and_destroy(sys->indepvars);  /* and indepvars */  
   }  
   if (sys->y != NULL && !sys->ycount) {  
     ascfree(sys->y);  
   }  
   if (sys->ydot != NULL && !sys->ydotcount) {  
     ascfree(sys->ydot);  
   }  
   if (sys->y_id != NULL) {  
     ascfree(sys->y_id);  
   }  
   if (sys->obs != NULL && !sys->obscount) {  
     ascfree(sys->obs);  
   }  
   if (sys->obs_id != NULL) {  
     ascfree(sys->obs_id);  
   }  
   ascfree(sys);  
 }  
   
 static  
 void IntegInitSymbols(void)  
 {  
   STATEFLAG = AddSymbol("ode_type");  
   STATEINDEX = AddSymbol("ode_id");  
   OBSINDEX = AddSymbol("obs_id");  
 }  
   
 /*  
  *  Returns a pointer to an ok struct Integ_system_t. Returns NULL  
  *  if one cannot be built from the instance given.  
  *  Diagnostics will be printed.  
  *  The pointer returned will also be stored in l_isys if someone  
  *  in Integrators.c needs it. Anyone who takes a copy of l_isys should  
  *  then set l_isys to NULL and keep track of the sys from there.  
  *  On exit the gl_list_t * in the returned pointer will all be NULL.  
  */  
 static struct Integ_system_t *Integ_system_build(CONST slv_system_t slvsys)  
 {  
   struct Integ_system_t *sys;  
   struct Integ_var_t *v1,*v2;  
   struct var_variable **vlist;  
   long half,i,len,vlen;  
   int happy=1;  
   
   if (slvsys==NULL) {  
     return NULL;  
   }  
   sys=(struct Integ_system_t *)asccalloc(1,sizeof(struct Integ_system_t));  
   if (sys==NULL) {  
     return sys;  
   }  
   l_isys = sys;  
   IntegInitSymbols();  
   
   /* collect potential states and derivatives */  
   sys->indepvars = gl_create(10L);  /* t var info */  
   sys->dynvars = gl_create(200L);  /* y ydot var info */  
   sys->obslist = gl_create(100L);  /* obs info */  
   if (sys->dynvars == NULL ||  
       sys->obslist == NULL ||  
       sys->indepvars == NULL ) {  
     FPRINTF(ASCERR,"ERROR: (Integ_system_build) Insufficient memory.\n");  
     Integ_system_destroy(sys);  
     l_isys = NULL;  
     return l_isys;  
   }  
   l_nstates = l_nderivs = 0;  
   /* visit all the slv_system_t master var lists to collect vars */  
   /* find the vars mostly in this one */  
   vlist = slv_get_master_var_list(slvsys);  
   vlen = slv_get_num_master_vars(slvsys);  
   for (i=0;i<vlen;i++) {  
     Integ_classify_vars(vlist[i]);  
   }  
   /* probably nothing here, but gotta check. */  
   vlist = slv_get_master_par_list(slvsys);  
   vlen = slv_get_num_master_pars(slvsys);  
   for (i=0;i<vlen;i++) {  
     Integ_classify_vars(vlist[i]);  
   }  
   /* might find t here */  
   vlist = slv_get_master_unattached_list(slvsys);  
   vlen = slv_get_num_master_unattached(slvsys);  
   for (i=0;i<vlen;i++) {  
     Integ_classify_vars(vlist[i]);  
   }  
   
   /* check the sanity of the independent variable */  
   len = gl_length(sys->indepvars);  
   if (!len) {  
     FPRINTF(ASCERR,  
       "ERROR: (Integ_system_build) No independent variable found.\n");  
     Integ_system_destroy(sys);  
     l_isys = NULL;  
     return l_isys;  
   }  
   if (len > 1) {  
     char *name;  
     FPRINTF(ASCERR,  
       "ERROR: (Integ_system_build) Excess %ld independent variables found:\n",  
       len);  
     for (i=1; i <=len;i++) {  
       v1 = (struct Integ_var_t *)gl_fetch(sys->indepvars,i);  
       if (v1 != NULL) {  
         name = WriteInstanceNameString(var_instance(v1->i),g_solvinst_cur);  
         if (name !=NULL) {  
           FPRINTF(ASCERR,"\t\n");  
           ascfree(name);  
           name = NULL;  
         }  
       }  
     }  
     FPRINTF(ASCERR, "Set the %s on all but one of these to %s >= 0.\n",  
             SCP(STATEFLAG),SCP(STATEFLAG));  
     Integ_system_destroy(sys);  
     l_isys = NULL;  
     return l_isys;  
   } else {  
     v1 = (struct Integ_var_t *)gl_fetch(sys->indepvars,1);  
     sys->x = v1->i;  
   }  
   /* check sanity of state and var lists */  
   len = gl_length(sys->dynvars);  
   half = len/2;  
   if (len % 2 || len == 0L || l_nstates != l_nderivs ) {  
     /* list length must be even for vars to pair off */  
     FPRINTF(ASCERR,"ERROR: (Integ_system_build) n_y != n_ydot\n");  
     FPRINTF(ASCERR,"\tor no dynamic vars found. Fix your indexing.\n");  
     Integ_system_destroy(sys);  
     l_isys = NULL;  
     return l_isys;  
   }  
   gl_sort(sys->dynvars,(CmpFunc)Integ_CmpDynVars);  
   if (gl_fetch(sys->dynvars,len)==NULL) {  
     FPRINTF(ASCERR,"ERROR: (Integ_system_build) Mysterious NULL found.\n");  
     FPRINTF(ASCERR,"\tPlease Report this to %s.\n",ASC_BIG_BUGMAIL);  
     Integ_system_destroy(sys);  
     l_isys = NULL;  
     return l_isys;  
   }  
   sys->states = gl_create(half);   /* state vars Integ_var_t references */  
   sys->derivs = gl_create(half);   /* derivative var atoms */  
   for (i=1;i < len; i+=2) {  
     v1 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i);  
     v2 = (struct Integ_var_t *)gl_fetch(sys->dynvars,i+1);  
     if (v1->type!=1  || v2 ->type !=2 || v1->index != v2->index) {  
       FPRINTF(ASCERR,"ERROR: (Integ_system_build) Mistyped or misindexed\n");  
       FPRINTF(ASCERR,  
              "\tdynamic variables: (%s = %ld,%s = %ld),(%s = %ld,%s = %ld).\n",  
              SCP(STATEFLAG),v1->type,SCP(STATEINDEX),v1->index,  
              SCP(STATEFLAG),v2->type,SCP(STATEINDEX),v2->index);  
       happy=0;  
       break;  
     } else {  
       gl_append_ptr(sys->states,(POINTER)v1);  
       gl_append_ptr(sys->derivs,(POINTER)v2->i);  
     }  
   }  
   if (!happy) {  
     Integ_system_destroy(sys);  
     l_isys = NULL;  
     return l_isys;  
   }  
   sys->n_y = half;  
   sys->y = (struct var_variable **)  
     ascmalloc(sizeof(struct var_variable *)*half);  
   sys->y_id = (long *)ascmalloc(sizeof(long)*half);  
   sys->ydot = (struct var_variable **)  
     ascmalloc(sizeof(struct var_variable *)*half);  
   if (sys->y==NULL || sys->ydot==NULL || sys->y_id==NULL) {  
     FPRINTF(ASCERR,"ERROR: (Integ_system_build) Insufficient memory.\n");  
     Integ_system_destroy(sys);  
     l_isys = NULL;  
     return l_isys;  
   }  
   for (i = 1; i <= half; i++) {  
     v1 = (struct Integ_var_t *)gl_fetch(sys->states,i);  
     sys->y[i-1] = v1->i;  
     sys->y_id[i-1] = v1->index;  
     sys->ydot[i-1] = (struct var_variable *)gl_fetch(sys->derivs,i);  
   }  
   
   /* reindex observations. sort if the user mostly numbered. take  
    * natural order if user just booleaned.  
    */  
   len = gl_length(sys->obslist);  
   /* we shouldn't be seeing NULL here ever except if malloc fail. */  
   if (len > 1L) {  
     half = ((struct Integ_var_t *)gl_fetch(sys->obslist,1))->index;  
     /* half != 0 now because we didn't collect 0 indexed vars */  
     for (i=2; i <= len; i++) {  
       if (half != ((struct Integ_var_t *)gl_fetch(sys->obslist,i))->index) {  
         /* change seen. sort and go on */  
         gl_sort(sys->obslist,(CmpFunc)Integ_CmpObs);  
         break;  
       }  
     }  
   }  
   for (i = half = 1; i <= len; i++) {  
     v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);  
     if (v2==NULL) {  
       /* we shouldn't be seeing NULL here ever except if malloc fail. */  
       gl_delete(sys->obslist,i,0); /* should not be gl_delete(so,i,1) */  
     } else {  
       Integ_SetObsId(v2->i,half);  
       v2->index = half++;  
     }  
   }  
   
   /* obslist now uniquely indexed, no nulls */  
   /* make into arrays */  
   half = gl_length(sys->obslist);  
   sys->obs = (struct var_variable **)  
     ascmalloc(sizeof(struct var_variable *)*half);  
   sys->obs_id = (long *)ascmalloc(sizeof(long)*half);  
   if ( sys->obs==NULL || sys->obs_id==NULL) {  
     FPRINTF(ASCERR,"ERROR: (Integ_system_build) Insufficient memory.\n");  
     Integ_system_destroy(sys);  
     l_isys = NULL;  
     return l_isys;  
   }  
   sys->n_obs = half;  
   for (i = 1; i <= half; i++) {  
     v2 = (struct Integ_var_t *)gl_fetch(sys->obslist,i);  
     sys->obs[i-1] = v2->i;  
     sys->obs_id[i-1] = v2->index;  
   }  
   
   /* don't need the gl_lists now that we have arrays for everyone */  
   gl_destroy(sys->states);  
   gl_destroy(sys->derivs);  
   gl_free_and_destroy(sys->indepvars);  /* we own the objects in indepvars */  
   gl_free_and_destroy(sys->dynvars);    /* we own the objects in dynvars */  
   gl_free_and_destroy(sys->obslist);    /* and obslist */  
   sys->states = NULL;  
   sys->derivs = NULL;  
   sys->indepvars = NULL;  
   sys->dynvars = NULL;  
   sys->obslist = NULL;  
   return l_isys;  
 }  
   
 /* boolean return. 0 = bad 1 = ok. BLSODE required setup stuff. */  
 static int Integ_setup_blsode(enum Integrator_type integrator)  
 {  
   struct Integ_system_t *blsys=NULL;  
   
   g_intgsys_cur = g_solvsys_cur;  
   if (g_intgsys_cur == NULL) {  
     FPRINTF(stderr, "g_intgsys_cur not correctly assigned.\n");  
     return 0;  
   }  
   /* verify integrator type ok. always passes for nonNULL inst. */  
   if (!(Asc_IntegInstIntegrable(g_solvinst_cur, integrator))) {  
     FPRINTF(stderr, "System is of wrong type to be ");  
     FPRINTF(stderr, "integrated using BLSODE.\n");  
     return 0;  
   }  
   
   /* this is a lie, but we will keep it.  
      We handle any linsol/linsolqr based solver. */  
   if (strcmp(slv_solver_name(slv_get_selected_solver(g_intgsys_cur)),"QRSlv")  
       != 0) {  
     FPRINTF(stderr, "QRSlv must be selected before integration.\n");  
     return 0;  
   }  
   
   blsys = Integ_system_build(g_solvsys_cur);  
   if (blsys == NULL) {  
     FPRINTF(ASCERR,"Unable to build blsode system.\n");  
     return 0;  
   }  
   
   g_intginst_d_n_eq = blsys->n_y;  
   g_intginst_d_n_obs = blsys->n_obs;  
   
   g_intginst_d_x = blsys->x;  
   g_intginst_d_y_vars = blsys->y;  
   g_intginst_d_dydx_vars = blsys->ydot;  
   g_intginst_d_obs_vars = blsys->obs;  
   blsys->ycount++;  
   blsys->ydotcount++;  
   blsys->obscount++;  
   
   print_debugstring("After values\n");  
   return 1;  
 }  
   
 /*  
  *  takes the type of integrator and sets up the global variables into  
  *  the current integration instance.  
  */  
 static int Integ_setup(enum Integrator_type integrator,  
                        long i0, long i1,  
                        double dt0,double dtmin,double dtmax,  
                        int moststeps)  
 {  
   long nstep;  
   unsigned long start_index=0, finish_index=0;  
   struct Integ_system_t *blsys;  
   
   switch (integrator) {  
   case BLSODE:  
     if (!Integ_setup_blsode(integrator)) {  
       return 0;  
     }  
     blsys = l_isys;  
     l_isys = NULL;  
     break;  
   default:  
     FPRINTF(stderr, "The requested type of integrator is ");  
     FPRINTF(stderr, "not currently available\n.");  
     return 0;  
   } /* END of integrator CASE */  
   
   switch (integrator) {  
   case BLSODE:  
     nstep = Asc_IntegGetNSamples()-1;  
     /* check for at least 2 steps and dimensionality of x vs steps here */  
     break;  
   default:  
     nstep = 0;  
     break;  
   }  
   if (i0<0 || i1 <0) {  
     FPRINTF(stdout, "An integration interval had been defined ");  
     FPRINTF(stdout, "from x[0] to x[%li].\n", nstep);  
     FPRINTF(stdout, "Enter start index: ");  
     start_index = (int)readlong((int)start_index);  
   
     if (start_index >= (unsigned long)nstep) {  
       FPRINTF(stderr, "ERROR: Start point < 0 OR start point >= %li.\n",nstep);  
       Integ_system_destroy(blsys);  
       return 0;  
     }  
   
     FPRINTF(stdout, "Enter finish index: ");  
     finish_index = readlong((int)finish_index);  
   
     if (finish_index > (unsigned long)nstep) {  
       FPRINTF(stderr, "ERROR: finish point < 0 OR finish point > %li.\n",nstep);  
       FPRINTF(stderr, "       finish point = %lu\n",finish_index);  
       Integ_system_destroy(blsys);  
       return 0;  
     }  
   } else {  
     start_index=i0;  
     finish_index =i1;  
     if (start_index >= (unsigned long)nstep) {  
       FPRINTF(stderr, "ERROR: Start point < 0 OR start point >= %li.\n",nstep);  
       FPRINTF(stderr, "       Start point = %lu\n",start_index);  
       Integ_system_destroy(blsys);  
       return 0;  
     }  
     if (finish_index > (unsigned long)nstep) {  
       FPRINTF(stderr, "ERROR: finish point < 0 OR finish point > %li.\n",nstep);  
       FPRINTF(stderr, "       finish point = %lu\n",finish_index);  
       Integ_system_destroy(blsys);  
       return 0;  
     }  
   }  
   if ((finish_index <= start_index) || (start_index >= finish_index)) {  
     FPRINTF(stderr, "ERROR: Finish point <= start point OR ");  
     FPRINTF(stderr, "start point >= finish point.\n");  
     FPRINTF(stderr, "       Start point = %lu\n",start_index);  
     FPRINTF(stderr, "       finish point = %lu\n",finish_index);  
     Integ_system_destroy(blsys);  
     return 0;  
   }  
   blsys->maxsteps = moststeps;  
   blsys->stepzero = dt0;  
   blsys->stepmin = dtmin;  
   blsys->stepmax = dtmax;  
   
   switch (integrator) {  
   case BLSODE:  
   /*  
    * Set up the Jacobian for the system. Does not check  
    * This assumes a clean sytem. Shut down will ensure a  
    * clean system at the end.  
    */  
     if (Integ_SetUpDiffs_BLsode(blsys)) {  
       return 0;  
     }  
     break;  
   default:  
     break;  
   }  
   Integ_Solve(integrator, start_index, finish_index,blsys);  
   Integ_system_destroy(blsys);  
   return 0;  
 }  
   
 /*  
  * switches:  
  * -engine $name  
  * -i0 $stepindex  
  * -i1 $stepindex  
  * -dt0 $initstepsize  
  * -dtmin $minstep  
  * -dtmax $maxstep  
  * -moststeps $moststeps  
  */  
 int Asc_IntegSetupCmd(ClientData cdata,Tcl_Interp *interp,  
                    int argc, CONST84 char *argv[])  
 {  
   enum Integrator_type integrator = UNKNOWN;  
   char buf[MAXIMUM_NUMERIC_LENGTH];         /* string to hold integer */  
   CONST84 char *engine = NULL;  
   int result = 0;         /* 0 = FALSE; 1 = TRUE */  
   long i0=(-1), i1=(-1);  
   int ifound = 0;  
   int k;  
   int moststeps=0;  
   double dt0=0, dtmin=0, dtmax=0;  
   CONST84 char *cdt0=NULL, *cdtmin=NULL, *cdtmax=NULL, *cmoststeps=NULL,  
        *ci0=NULL, *ci1=NULL;  
   
   (void)cdata;    /* stop gcc whine about unused parameter */  
   
   k = 1;  
   while (k < (argc-1)) { /* arguments come in pairs */  
     if (strcmp(argv[k],"-engine")==0) {  
       engine = argv[k+1];  
       k+=2;  
       continue;  
     }  
     if (strcmp(argv[k],"-i1")==0) {  
       ci1 = argv[k+1];  
       k+=2;  
       continue;  
     }  
     if (strcmp(argv[k],"-i0")==0) {  
       ci0 = argv[k+1];  
       k+=2;  
       continue;  
     }  
     if (strcmp(argv[k],"-moststeps")==0) {  
       cmoststeps = argv[k+1];  
       k+=2;  
       continue;  
     }  
     if (strcmp(argv[k],"-dtmax")==0) {  
       cdtmax = argv[k+1];  
       k+=2;  
       continue;  
     }  
     if (strcmp(argv[k],"-dtmin")==0) {  
       cdtmin = argv[k+1];  
       k+=2;  
       continue;  
     }  
     if (strcmp(argv[k],"-dt0")==0) {  
       cdt0 = argv[k+1];  
       k+=2;  
       continue;  
     }  
     Tcl_AppendResult(interp,argv[0],": unrecognized option ",  
                      argv[k],".",SNULL);  
     return TCL_ERROR;  
   }  
   
   if (engine != NULL && strncmp(engine,"BLSODE",3)==0) {  
     integrator = BLSODE;  
     ifound=1;  
   }  
   if (!ifound) {  
     Tcl_SetResult(interp, "Unsupported integrator", TCL_STATIC);  
     Tcl_AppendResult(interp," ",engine,SNULL);  
     return TCL_ERROR;  
   }  
   if (ci0 != NULL && ci1 != NULL) {  
     /* get i0, i1 if both supplied. */  
     long i;  
     if (Tcl_ExprLong(interp,ci0,&i)==TCL_ERROR|| i<0) {  
       Tcl_ResetResult(interp);  
       Tcl_SetResult(interp, "integrate_setup: index i0 invalid", TCL_STATIC);  
       return TCL_ERROR;  
     }  
     i0=i;  
     if (Tcl_ExprLong(interp,ci1,&i)==TCL_ERROR|| i<i0) {  
       Tcl_ResetResult(interp);  
       Tcl_SetResult(interp, "integrate_setup: index i1 invalid", TCL_STATIC);  
       return TCL_ERROR;  
     }  
     i1=i;  
   }  
   if (cdt0 != NULL) {  
     if (Tcl_GetDouble(interp,cdt0,&dt0) != TCL_OK) {  
       Tcl_ResetResult(interp);  
       Tcl_AppendResult(interp, "integrate_setup: initial step length invalid",  
                        " (",cdt0,")", SNULL);  
       return TCL_ERROR;  
     }  
   }  
   if (cdtmin != NULL) {  
     if (Tcl_GetDouble(interp,cdtmin,&dtmin) != TCL_OK || dtmin < 0) {  
       Tcl_ResetResult(interp);  
       Tcl_AppendResult(interp, "integrate_setup: minimum step length invalid",  
                        " (",cdtmin,")", SNULL);  
       return TCL_ERROR;  
     }  
   }  
   if (cdtmax != NULL) {  
     if (Tcl_GetDouble(interp,cdtmax,&dtmax) != TCL_OK || dtmax < dtmin) {  
       Tcl_ResetResult(interp);  
       Tcl_AppendResult(interp, "integrate_setup: maximum step length invalid",  
                        " (",cdtmax,")", SNULL);  
       return TCL_ERROR;  
     }  
   }  
   if (cmoststeps != NULL) {  
     if (Tcl_GetInt(interp,cmoststeps,&moststeps) != TCL_OK || moststeps < 0) {  
       Tcl_ResetResult(interp);  
       Tcl_AppendResult(interp, "integrate_setup: maximum internal steps bad",  
                        " (",cmoststeps,")", SNULL);  
       return TCL_ERROR;  
     }  
   }  
   result = Integ_setup(integrator,i0,i1,dt0,dtmin,dtmax,moststeps);  
   sprintf(buf, "%d", result);  
   Tcl_SetResult(interp, buf, TCL_VOLATILE);  
   return TCL_OK;  
 }  
   
 /*  
  *  Deallocates any memory used and sets all integration global points  
  *  to NULL.  
  */  
 static int Integ_Cleanup() {  
   g_intgsys_cur = NULL;  
   g_intginst_d_x = NULL;  
   if(g_intginst_d_y_vars !=NULL) {  
     ascfree(g_intginst_d_y_vars);  
   }  
   if(g_intginst_d_dydx_vars !=NULL) {  
     ascfree(g_intginst_d_dydx_vars);  
   }  
   if(g_intginst_d_obs_vars !=NULL) {  
     ascfree(g_intginst_d_obs_vars);  
   }  
   g_intginst_d_y_vars=NULL;  
   g_intginst_d_dydx_vars=NULL;  
   g_intginst_d_obs_vars=NULL;  
   
   /*  
    * Cleanup the derivative stuff.  
    */  
   if (g_intg_diff.input_indices) {  
     ascfree((char *)g_intg_diff.input_indices);  
   }  
   if (g_intg_diff.output_indices) {  
     ascfree((char *)g_intg_diff.output_indices);  
   }  
   if (g_intg_diff.y_vars) {  
     ascfree((char *)g_intg_diff.y_vars);  
   }  
   if (g_intg_diff.ydot_vars) {  
     ascfree((char *)g_intg_diff.ydot_vars);  
   }  
   if (g_intg_diff.rlist) {  
     ascfree((char *)g_intg_diff.rlist);  
   }  
   if (g_intg_diff.dydx_dx) {  
     DestroyDenseMatrix(g_intg_diff.dydx_dx, g_intg_diff.n_eqns);  
   }  
   
   g_intg_diff.input_indices = NULL;  
   g_intg_diff.output_indices = NULL;  
   g_intg_diff.y_vars = NULL;  
   g_intg_diff.ydot_vars = NULL;  
   g_intg_diff.dydx_dx =  NULL;  
   g_intg_diff.rlist =  NULL;  
   g_intg_diff.n_eqns = 0L;  
   
   print_debugstring("integrate_cleanup\n");  
   return 0;  
 }  
   
 /********************************************************************/  
   
 int Asc_IntegCleanupCmd(ClientData cdata,Tcl_Interp *interp,  
                      int argc, CONST84 char *argv[])  
 {  
   (void)cdata;    /* stop gcc whine about unused parameter */  
   (void)argv;     /* stop gcc whine about unused parameter */  
   
   if (argc!=1) {  
     Tcl_SetResult(interp, "integrate_cleanup takes no arguments", TCL_STATIC);  
     return TCL_ERROR;  
   }  
   
   Integ_Cleanup();  
   return TCL_OK;  
 }  
Legend:
Removed from v.571  
changed lines
  Added in v.670
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22