/[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 689 by johnpye, Wed Jun 21 07:57:19 2006 UTC revision 690 by johnpye, Thu Jun 22 00:48:31 2006 UTC
# Line 1  Line 1 
1  /*  ASCEND modelling environment  /*  ASCEND modelling environment
2      Copyright 1997, Carnegie Mellon University      Copyright 1997, Carnegie Mellon University
3      Copyright (C) 2006 Carnegie Mellon University      Copyright (C) 2006 Carnegie Mellon University
4    
5      This program is free software; you can redistribute it and/or modify      This program is free software; you can redistribute it and/or modify
6      it under the terms of the GNU General Public License as published by      it under the terms of the GNU General Public License as published by
7      the Free Software Foundation; either version 2, or (at your option)      the Free Software Foundation; either version 2, or (at your option)
8      any later version.      any later version.
9    
10      This program is distributed in the hope that it will be useful,      This program is distributed in the hope that it will be useful,
11      but WITHOUT ANY WARRANTY; without even the implied warranty of      but WITHOUT ANY WARRANTY; without even the implied warranty of
12      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13      GNU General Public License for more details.      GNU General Public License for more details.
14    
15      You should have received a copy of the GNU General Public License      You should have received a copy of the GNU General Public License
16      along with this program; if not, write to the Free Software      along with this program; if not, write to the Free Software
17      Foundation, Inc., 59 Temple Place - Suite 330,      Foundation, Inc., 59 Temple Place - Suite 330,
18      Boston, MA 02111-1307, USA.      Boston, MA 02111-1307, USA.
19  *//**  *//**
20      @file      @file
21      Tcl/Tk interface functions for the Integration feature      Tcl/Tk interface functions for the Integration feature
22  *//*  *//*
23      by Kirk Abbott and Ben Allan      by Kirk Abbott and Ben Allan
24      Created: 1/94      Created: 1/94
25      Version: $Revision: 1.32 $      Version: $Revision: 1.32 $
26      Version control file: $RCSfile: Integrators.c,v $      Version control file: $RCSfile: Integrators.c,v $
27      Date last modified: $Date: 2003/08/23 18:43:06 $      Date last modified: $Date: 2003/08/23 18:43:06 $
28      Last modified by: $Author: ballan $      Last modified by: $Author: ballan $
29  */  */
30    
31  #include <tcl.h>  #include <tcl.h>
32    
33  #include <utilities/ascConfig.h>  #include <utilities/ascConfig.h>
34  #include <solver/integrator.h>  #include <solver/integrator.h>
35  #include <solver/lsode.h>  #include <solver/lsode.h>
36  #include <solver/samplelist.h>  #include <solver/samplelist.h>
37    
38  #include "HelpProc.h"  #include "HelpProc.h"
39  #include "Integrators.h"  #include "Integrators.h"
40  #include "BrowserQuery.h"  #include "BrowserQuery.h"
41  #include "Qlfdid.h"  #include "Qlfdid.h"
42  #include "UnitsProc.h"  #include "UnitsProc.h"
43  #include "BrowserProc.h"  #include "BrowserProc.h"
44  #include "HelpProc.h"  #include "HelpProc.h"
45  #include "SolverGlobals.h"  #include "SolverGlobals.h"
46    
47  #ifndef lint  #ifndef lint
48  static CONST char IntegratorsID[] = "$Id: Integrators.c,v 1.32 2003/08/23 18:43:06 ballan Exp $";  static CONST char IntegratorsID[] = "$Id: Integrators.c,v 1.32 2003/08/23 18:43:06 ballan Exp $";
49  #endif  #endif
50    
51  #define SNULL (char *)NULL  #define SNULL (char *)NULL
52    
53  static SampleList l_samplelist;  static SampleList l_samplelist;
54    
55  /*-------------------------------------------------------  /*-------------------------------------------------------
56    HANDLING OF OUTPUT FILES    HANDLING OF OUTPUT FILES
57  */  */
58    
59  /* vars relating to output files */  /* vars relating to output files */
60  static FILE *l_obs_file = NULL;  static FILE *l_obs_file = NULL;
61  static char *l_obs_filename = NULL;  static char *l_obs_filename = NULL;
62    
63  static FILE *l_y_file = NULL;  static FILE *l_y_file = NULL;
64  static char *l_y_filename = NULL;  static char *l_y_filename = NULL;
65    
66  static int l_print_option = 1; /* default si */  static int l_print_option = 1; /* default si */
67  static int l_print_fmt = 0; /* default variable */  static int l_print_fmt = 0; /* default variable */
68    
69  FILE *Asc_IntegOpenYFile(void)  FILE *Asc_IntegOpenYFile(void)
70  {  {
71    if (l_y_filename==NULL) {    if (l_y_filename==NULL) {
72      return NULL;      return NULL;
73    }    }
74    l_y_file = fopen(l_y_filename,"a+");    l_y_file = fopen(l_y_filename,"a+");
75    
76    if (l_y_file==NULL) {    if (l_y_file==NULL) {
77      FPRINTF(ASCERR,      FPRINTF(ASCERR,
78        "WARNING: (integrate) Unable to open\n\t%s\nfor state output log.\n",        "WARNING: (integrate) Unable to open\n\t%s\nfor state output log.\n",
79        l_y_filename);        l_y_filename);
80    } else {    } else {
81      time_t t;      time_t t;
82      t = time((time_t *)NULL);      t = time((time_t *)NULL);
83      FPRINTF(l_y_file,"DATASET %s", asctime(localtime(&t)));      FPRINTF(l_y_file,"DATASET %s", asctime(localtime(&t)));
84      FFLUSH(l_y_file);      FFLUSH(l_y_file);
85    }    }
86    return l_y_file;    return l_y_file;
87  }  }
88    
89  FILE *Asc_IntegOpenObsFile(void)  FILE *Asc_IntegOpenObsFile(void)
90  {  {
91    if (l_obs_filename==NULL) {    if (l_obs_filename==NULL) {
92      return NULL;      return NULL;
93    }    }
94    l_obs_file = fopen(l_obs_filename,"a+");    l_obs_file = fopen(l_obs_filename,"a+");
95    if (l_obs_file==NULL) {    if (l_obs_file==NULL) {
96      FPRINTF(ASCERR,      FPRINTF(ASCERR,
97        "WARNING: (integrate) Unable to open\n\t%s\nfor observation log.\n",        "WARNING: (integrate) Unable to open\n\t%s\nfor observation log.\n",
98        l_obs_filename);        l_obs_filename);
99    } else {    } else {
100      time_t t;      time_t t;
101      t = time((time_t *)NULL);      t = time((time_t *)NULL);
102      FPRINTF(l_obs_file,"DATASET %s", asctime(localtime(&t)));      FPRINTF(l_obs_file,"DATASET %s", asctime(localtime(&t)));
103      FFLUSH(l_obs_file);      FFLUSH(l_obs_file);
104    }    }
105    return l_obs_file;    return l_obs_file;
106  }  }
107    
108  FILE *Asc_IntegGetYFile(void)  FILE *Asc_IntegGetYFile(void)
109  {  {
110    return l_y_file;    return l_y_file;
111  }  }
112    
113  FILE *Asc_IntegGetObsFile(void)  FILE *Asc_IntegGetObsFile(void)
114  {  {
115    return l_obs_file;    return l_obs_file;
116  }  }
117    
118  void Asc_IntegReleaseYFile(void)  void Asc_IntegReleaseYFile(void)
119  {  {
120    l_y_file = NULL;    l_y_file = NULL;
121  }  }
122  void Asc_IntegReleaseObsFile(void)  void Asc_IntegReleaseObsFile(void)
123  {  {
124    l_obs_file = NULL;    l_obs_file = NULL;
125  }  }
126    
127  /********************************************************************/  /********************************************************************/
128  /* string column layout: 26 char columns  if l_print_fmt else variable */  /* string column layout: 26 char columns  if l_print_fmt else variable */
129  /* numeric format: leading space plus characters from Unit*Value */  /* numeric format: leading space plus characters from Unit*Value */
130  /* index format: left padded long int */  /* index format: left padded long int */
131  /* headline format space plus - s */  /* headline format space plus - s */
132  #define BCOLSFMT (l_print_fmt ? "%-26s" : "\t%s")  #define BCOLSFMT (l_print_fmt ? "%-26s" : "\t%s")
133  #define BCOLNFMT (l_print_fmt ? " %-25s" : "\t%s")  #define BCOLNFMT (l_print_fmt ? " %-25s" : "\t%s")
134  #define BCOLIFMT (l_print_fmt ? " %25ld" : "\t%ld")  #define BCOLIFMT (l_print_fmt ? " %25ld" : "\t%ld")
135  #define BCOLHFMT (l_print_fmt ? " -------------------------" : "\t---")  #define BCOLHFMT (l_print_fmt ? " -------------------------" : "\t---")
136    
137  void Asc_IntegPrintYHeader(FILE *fp, IntegratorSystem *blsys)  void Asc_IntegPrintYHeader(FILE *fp, IntegratorSystem *blsys)
138  {  {
139    long i,len, *yip;    long i,len, *yip;
140    char *name;    char *name;
141    struct Instance *in;    struct Instance *in;
142    int si;    int si;
143    if (fp==NULL) {    if (fp==NULL) {
144      return;      return;
145    }    }
146    if (blsys==NULL) {    if (blsys==NULL) {
147      FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintYHeader: called w/o data\n");      FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintYHeader: called w/o data\n");
148      return;      return;
149    }    }
150    if (blsys->n_y == 0) {    if (blsys->n_y == 0) {
151      return;      return;
152    }    }
153    if (blsys->y == NULL) {    if (blsys->y == NULL) {
154      FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintYHeader: called w/NULL data\n");      FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintYHeader: called w/NULL data\n");
155      return;      return;
156    }    }
157    len = blsys->n_y;    len = blsys->n_y;
158    yip = blsys->y_id;    yip = blsys->y_id;
159    si = l_print_option;    si = l_print_option;
160    /* output indep var name */    /* output indep var name */
161    /* output dep var names */    /* output dep var names */
162    FPRINTF(fp,"States: (user index) (name) (units)\n");    FPRINTF(fp,"States: (user index) (name) (units)\n");
163    in = var_instance(blsys->x);    in = var_instance(blsys->x);
164    FPRINTF(fp,"{indvar}");  /* indep id name */    FPRINTF(fp,"{indvar}");  /* indep id name */
165    name = WriteInstanceNameString(in,g_solvinst_cur);    name = WriteInstanceNameString(in,g_solvinst_cur);
166    FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));    FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));
167    ascfree(name);    ascfree(name);
168    for (i=0; i< len; i++) {    for (i=0; i< len; i++) {
169      in = var_instance(blsys->y[i]);      in = var_instance(blsys->y[i]);
170      FPRINTF(fp,"{%ld}",yip[i]);  /* user id # */      FPRINTF(fp,"{%ld}",yip[i]);  /* user id # */
171      name = WriteInstanceNameString(in,g_solvinst_cur);      name = WriteInstanceNameString(in,g_solvinst_cur);
172      FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));      FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));
173      ascfree(name);      ascfree(name);
174    }    }
175    FPRINTF(fp,BCOLSFMT,"indvar");    FPRINTF(fp,BCOLSFMT,"indvar");
176    for (i=0; i < len; i++) {    for (i=0; i < len; i++) {
177      FPRINTF(fp,BCOLIFMT,yip[i]);      FPRINTF(fp,BCOLIFMT,yip[i]);
178    }    }
179    FPRINTF(fp,"\n");    FPRINTF(fp,"\n");
180    for (i=0; i <= len; i++) {    for (i=0; i <= len; i++) {
181      FPRINTF(fp,BCOLHFMT);      FPRINTF(fp,BCOLHFMT);
182    }    }
183    FPRINTF(fp,"\n");    FPRINTF(fp,"\n");
184  }  }
185  /********************************************************************/  /********************************************************************/
186  void Asc_IntegPrintObsHeader(FILE *fp, IntegratorSystem *blsys)  void Asc_IntegPrintObsHeader(FILE *fp, IntegratorSystem *blsys)
187  {  {
188    long i,len, *obsip;    long i,len, *obsip;
189    char *name;    char *name;
190    struct Instance *in;    struct Instance *in;
191    int si;    int si;
192    if (fp==NULL) {    if (fp==NULL) {
193      return;      return;
194    }    }
195    if (blsys==NULL) {    if (blsys==NULL) {
196      ERROR_REPORTER_HERE(ASC_PROG_ERR,"called without data");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"called without data");
197      return;      return;
198    }    }
199    if (blsys->n_obs == 0) {    if (blsys->n_obs == 0) {
200      return;      return;
201    }    }
202    if (blsys->obs == NULL) {    if (blsys->obs == NULL) {
203      ERROR_REPORTER_HERE(ASC_PROG_ERR,"called with NULL data");      ERROR_REPORTER_HERE(ASC_PROG_ERR,"called with NULL data");
204      return;      return;
205    }    }
206    len = blsys->n_obs;    len = blsys->n_obs;
207    obsip = blsys->obs_id;    obsip = blsys->obs_id;
208    si = l_print_option;    si = l_print_option;
209    FPRINTF(fp,"Observations: (user index) (name) (units)\n");    FPRINTF(fp,"Observations: (user index) (name) (units)\n");
210    /* output indep var name */    /* output indep var name */
211    /* output obs var names */    /* output obs var names */
212    in = var_instance(blsys->x);    in = var_instance(blsys->x);
213    FPRINTF(fp,"{indvar}");  /* indep id name */    FPRINTF(fp,"{indvar}");  /* indep id name */
214    name = WriteInstanceNameString(in,g_solvinst_cur);    name = WriteInstanceNameString(in,g_solvinst_cur);
215    FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));    FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));
216    ascfree(name);    ascfree(name);
217    for (i=0; i< len; i++) {    for (i=0; i< len; i++) {
218      in = var_instance(blsys->obs[i]);      in = var_instance(blsys->obs[i]);
219      FPRINTF(fp,"{%ld}",obsip[i]);  /* user id # */      FPRINTF(fp,"{%ld}",obsip[i]);  /* user id # */
220      name = WriteInstanceNameString(in,g_solvinst_cur);      name = WriteInstanceNameString(in,g_solvinst_cur);
221      FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));      FPRINTF(fp,"\t{%s}\t{%s}\n",name,Asc_UnitString(in,si));
222      ascfree(name);      ascfree(name);
223    }    }
224    FPRINTF(fp,BCOLSFMT,"indvar");    FPRINTF(fp,BCOLSFMT,"indvar");
225    for (i=0; i < len; i++) {    for (i=0; i < len; i++) {
226      FPRINTF(fp,BCOLIFMT,obsip[i]);      FPRINTF(fp,BCOLIFMT,obsip[i]);
227    }    }
228    FPRINTF(fp,"\n");    FPRINTF(fp,"\n");
229    for (i=0; i <= len; i++) {    for (i=0; i <= len; i++) {
230      FPRINTF(fp,BCOLHFMT);      FPRINTF(fp,BCOLHFMT);
231    }    }
232    FPRINTF(fp,"\n");    FPRINTF(fp,"\n");
233  }  }
234    
235  /********************************************************************/  /********************************************************************/
236  void Asc_IntegPrintYLine(FILE *fp, IntegratorSystem *blsys)  void Asc_IntegPrintYLine(FILE *fp, IntegratorSystem *blsys)
237  {  {
238    long i,len;    long i,len;
239    struct var_variable **vp;    struct var_variable **vp;
240    int si;    int si;
241    if (fp==NULL) {    if (fp==NULL) {
242      return;      return;
243    }    }
244    if (blsys==NULL) {    if (blsys==NULL) {
245      FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintYLine: called w/o data\n");      FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintYLine: called w/o data\n");
246      return;      return;
247    }    }
248    if (blsys->n_y == 0) {    if (blsys->n_y == 0) {
249      return;      return;
250    }    }
251    if (blsys->y == NULL) {    if (blsys->y == NULL) {
252      FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintYHeader: called w/NULL data\n");      FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintYHeader: called w/NULL data\n");
253      return;      return;
254    }    }
255    vp = blsys->y;    vp = blsys->y;
256    len = blsys->n_y;    len = blsys->n_y;
257    si = l_print_option;    si = l_print_option;
258    FPRINTF(fp,BCOLNFMT,Asc_UnitlessValue(var_instance(blsys->x),si));    FPRINTF(fp,BCOLNFMT,Asc_UnitlessValue(var_instance(blsys->x),si));
259    for (i=0; i < len; i++) {    for (i=0; i < len; i++) {
260      FPRINTF(fp,BCOLNFMT, Asc_UnitlessValue(var_instance(vp[i]),si));      FPRINTF(fp,BCOLNFMT, Asc_UnitlessValue(var_instance(vp[i]),si));
261    }    }
262    FPRINTF(fp,"\n");    FPRINTF(fp,"\n");
263  }  }
264    
265  void Asc_IntegPrintObsLine(FILE *fp, IntegratorSystem *blsys)  void Asc_IntegPrintObsLine(FILE *fp, IntegratorSystem *blsys)
266  {  {
267    long i,len;    long i,len;
268    struct var_variable **vp;    struct var_variable **vp;
269    int si;    int si;
270    if (fp==NULL) {    if (fp==NULL) {
271      return;      return;
272    }    }
273    if (blsys==NULL) {    if (blsys==NULL) {
274      FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintObsLine: called w/o data\n");      FPRINTF(ASCERR,"WARNING: (Asc_IntegPrintObsLine: called w/o data\n");
275      return;      return;
276    }    }
277    if (blsys->n_obs == 0) {    if (blsys->n_obs == 0) {
278      return;      return;
279    }    }
280    if (blsys->obs == NULL) {    if (blsys->obs == NULL) {
281      FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintObsHeader: called w/NULL data\n");      FPRINTF(ASCERR,"ERROR: (Asc_IntegPrintObsHeader: called w/NULL data\n");
282      return;      return;
283    }    }
284    vp = blsys->obs;    vp = blsys->obs;
285    len = blsys->n_obs;    len = blsys->n_obs;
286    si = l_print_option;    si = l_print_option;
287    FPRINTF(fp,BCOLNFMT,Asc_UnitlessValue(var_instance(blsys->x),si));    FPRINTF(fp,BCOLNFMT,Asc_UnitlessValue(var_instance(blsys->x),si));
288    for (i=0; i < len; i++) {    for (i=0; i < len; i++) {
289      FPRINTF(fp,BCOLNFMT, Asc_UnitlessValue(var_instance(vp[i]),si));      FPRINTF(fp,BCOLNFMT, Asc_UnitlessValue(var_instance(vp[i]),si));
290    }    }
291    FPRINTF(fp,"\n");    FPRINTF(fp,"\n");
292  }  }
293    
294    
295  /*---------------------------------------------*/  /*---------------------------------------------*/
296    
297  int Asc_IntegSetYFileCmd(ClientData cdata, Tcl_Interp *interp,  int Asc_IntegSetYFileCmd(ClientData cdata, Tcl_Interp *interp,
298                        int argc, CONST84 char *argv[])                        int argc, CONST84 char *argv[])
299  {  {
300    size_t len;    size_t len;
301    
302    UNUSED_PARAMETER(cdata);    UNUSED_PARAMETER(cdata);
303    
304    if ( argc != 2 ) {    if ( argc != 2 ) {
305      FPRINTF(ASCERR, "integrate_set_y_file: called without filename.\n");      FPRINTF(ASCERR, "integrate_set_y_file: called without filename.\n");
306      Tcl_SetResult(interp,      Tcl_SetResult(interp,
307                    "integrate_set_y_file <filename,""> called without arg.",                    "integrate_set_y_file <filename,""> called without arg.",
308                    TCL_STATIC);                    TCL_STATIC);
309      return TCL_ERROR;      return TCL_ERROR;
310    }    }
311    if (len >0 ) {    if (len >0 ) {
312      l_y_filename = Asc_MakeInitString((int)len);      l_y_filename = Asc_MakeInitString((int)len);
313      sprintf(l_y_filename,"%s",argv[1]);      sprintf(l_y_filename,"%s",argv[1]);
314    } else {    } else {
315      l_y_filename = NULL;      l_y_filename = NULL;
316    }    }
317    return TCL_OK;    return TCL_OK;
318  }  }
319    
320  /*---------------------------------------------*/  /*---------------------------------------------*/
321    
322  int Asc_IntegSetObsFileCmd(ClientData cdata, Tcl_Interp *interp,  int Asc_IntegSetObsFileCmd(ClientData cdata, Tcl_Interp *interp,
323                          int argc, CONST84 char *argv[])                          int argc, CONST84 char *argv[])
324  {  {
325    size_t len;    size_t len;
326    
327    UNUSED_PARAMETER(cdata);    UNUSED_PARAMETER(cdata);
328    
329    if ( argc != 2 ) {    if ( argc != 2 ) {
330      FPRINTF(ASCERR, "integrate_set_obs_file: called without filename.\n");      FPRINTF(ASCERR, "integrate_set_obs_file: called without filename.\n");
331      Tcl_SetResult(interp,      Tcl_SetResult(interp,
332                    "integrate_set_obs_file <filename,""> called without arg.",                    "integrate_set_obs_file <filename,""> called without arg.",
333                    TCL_STATIC);                    TCL_STATIC);
334      return TCL_ERROR;      return TCL_ERROR;
335    }    }
336    if (l_obs_filename != NULL) {    if (l_obs_filename != NULL) {
337      ascfree(l_obs_filename);      ascfree(l_obs_filename);
338    }    }
339    len = strlen(argv[1]);    len = strlen(argv[1]);
340    if (len >0 ) {    if (len >0 ) {
341      l_obs_filename = Asc_MakeInitString((int)len);      l_obs_filename = Asc_MakeInitString((int)len);
342      sprintf(l_obs_filename,"%s",argv[1]);      sprintf(l_obs_filename,"%s",argv[1]);
343    } else {    } else {
344      l_obs_filename = NULL;      l_obs_filename = NULL;
345    }    }
346    return TCL_OK;    return TCL_OK;
347  }  }
348    
349  /*---------------------------------------------*/  /*---------------------------------------------*/
350    
351  int Asc_IntegSetFileUnitsCmd(ClientData cdata, Tcl_Interp *interp,  int Asc_IntegSetFileUnitsCmd(ClientData cdata, Tcl_Interp *interp,
352                              int argc, CONST84 char *argv[])                              int argc, CONST84 char *argv[])
353  {  {
354    UNUSED_PARAMETER(cdata);    UNUSED_PARAMETER(cdata);
355    
356    if ( argc != 2 ) {    if ( argc != 2 ) {
357      FPRINTF(ASCERR, "integrate_logunits: called without printoption.\n");      FPRINTF(ASCERR, "integrate_logunits: called without printoption.\n");
358      Tcl_SetResult(interp,"integrate_logunits <display,si> called without arg.",      Tcl_SetResult(interp,"integrate_logunits <display,si> called without arg.",
359                    TCL_STATIC);                    TCL_STATIC);
360      return TCL_ERROR;      return TCL_ERROR;
361    }    }
362    switch (argv[1][0]) {    switch (argv[1][0]) {
363    case 's':    case 's':
364      l_print_option = 1;      l_print_option = 1;
365      break;      break;
366    case 'd':    case 'd':
367      l_print_option = 0;      l_print_option = 0;
368      break;      break;
369    default:    default:
370      FPRINTF(ASCERR,"integrate_logunits: called with bogus argument.\n");      FPRINTF(ASCERR,"integrate_logunits: called with bogus argument.\n");
371      FPRINTF(ASCERR,"logunits remain set to %s.\n",      FPRINTF(ASCERR,"logunits remain set to %s.\n",
372        (l_print_option ? "si":"display"));        (l_print_option ? "si":"display"));
373      break;      break;
374    }    }
375    return TCL_OK;    return TCL_OK;
376  }  }
377    
378  /*---------------------------------------------*/  /*---------------------------------------------*/
379    
380  int Asc_IntegSetFileFormatCmd(ClientData cdata, Tcl_Interp *interp,  int Asc_IntegSetFileFormatCmd(ClientData cdata, Tcl_Interp *interp,
381                               int argc, CONST84 char *argv[])                               int argc, CONST84 char *argv[])
382  {  {
383    UNUSED_PARAMETER(cdata);    UNUSED_PARAMETER(cdata);
384    
385    if ( argc != 2 ) {    if ( argc != 2 ) {
386      FPRINTF(ASCERR, "integrate_logformat called without printoption.\n");      FPRINTF(ASCERR, "integrate_logformat called without printoption.\n");
387      Tcl_SetResult(interp,      Tcl_SetResult(interp,
388                    "integrate_logformat <fixed,variable> called without arg.",                    "integrate_logformat <fixed,variable> called without arg.",
389                    TCL_STATIC);                    TCL_STATIC);
390      return TCL_ERROR;      return TCL_ERROR;
391    }    }
392    switch (argv[1][0]) {    switch (argv[1][0]) {
393    case 'f':    case 'f':
394      l_print_fmt = 1;      l_print_fmt = 1;
395      break;      break;
396    case 'v':    case 'v':
397      l_print_fmt = 0;      l_print_fmt = 0;
398      break;      break;
399    default:    default:
400      FPRINTF(ASCERR,"integrate_logformat: called with bogus argument.\n");      FPRINTF(ASCERR,"integrate_logformat: called with bogus argument.\n");
401      FPRINTF(ASCERR,"logformat remains set to %s.\n",      FPRINTF(ASCERR,"logformat remains set to %s.\n",
402        (l_print_fmt ? "fixed":"variable"));        (l_print_fmt ? "fixed":"variable"));
403      break;      break;
404    }    }
405    return TCL_OK;    return TCL_OK;
406  }  }
407    
408    
409  /*---------------------------------------------------------------  /*---------------------------------------------------------------
410    DEFINING THE TIMESTEPS FOR INTEGRATION    DEFINING THE TIMESTEPS FOR INTEGRATION
411  */  */
412    
413  int Asc_IntegGetXSamplesCmd(ClientData cdata, Tcl_Interp *interp,  int Asc_IntegGetXSamplesCmd(ClientData cdata, Tcl_Interp *interp,
414                            int argc, CONST84 char *argv[])                            int argc, CONST84 char *argv[])
415  {  {
416    static char sval[40]; /* buffer long enough to hold a double printed */    static char sval[40]; /* buffer long enough to hold a double printed */
417    struct Units *du = NULL;    struct Units *du = NULL;
418    const dim_type *dp;    const dim_type *dp;
419    long i,len;    long i,len;
420    double *uvalues = NULL;    double *uvalues = NULL;
421    char *ustring;    char *ustring;
422    double *uv;    double *uv;
423    int trydu=0, prec, stat=0;    int trydu=0, prec, stat=0;
424    SampleList *samples;    SampleList *samples;
425    
426    UNUSED_PARAMETER(cdata);    UNUSED_PARAMETER(cdata);
427    
428    if (( argc < 1 ) || ( argc > 2 )) {    if (( argc < 1 ) || ( argc > 2 )) {
429      Tcl_SetResult(interp,      Tcl_SetResult(interp,
430                    "integrate_get_samples: expected 0 or 1 args [display]",                    "integrate_get_samples: expected 0 or 1 args [display]",
431                    TCL_STATIC);                    TCL_STATIC);
432      return TCL_ERROR;      return TCL_ERROR;
433    }    }
434    if (argc==2) {    if (argc==2) {
435      trydu = 1;      trydu = 1;
436      if( argv[1][0] != 'd') {      if( argv[1][0] != 'd') {
437        Tcl_SetResult(interp, "integrate_get_samples: expected display but got ",        Tcl_SetResult(interp, "integrate_get_samples: expected display but got ",
438                      TCL_STATIC);                      TCL_STATIC);
439        Tcl_AppendResult(interp,argv[1],".",SNULL);        Tcl_AppendResult(interp,argv[1],".",SNULL);
440        return TCL_ERROR;        return TCL_ERROR;
441      }      }
442    }    }
443    
444    len = samplelist_length(&l_samplelist);    len = samplelist_length(&l_samplelist);
445    dp = samplelist_dim(&l_samplelist);    dp = samplelist_dim(&l_samplelist);
446    
447    if (len <1) {    if (len <1) {
448      Tcl_SetResult(interp, "{} {}", TCL_STATIC);      Tcl_SetResult(interp, "{} {}", TCL_STATIC);
449      return TCL_OK;      return TCL_OK;
450    }    }
451    
452    if (trydu) {    if (trydu) {
453    
454      /* Allocate the space for the retrieved values... */      /* Allocate the space for the retrieved values... */
455      uvalues = ASC_NEW_ARRAY(double,len);      uvalues = ASC_NEW_ARRAY(double,len);
456      if (uvalues == NULL) {      if (uvalues == NULL) {
457        Tcl_SetResult(interp, "integrate_get_samples: Insufficient memory.",        Tcl_SetResult(interp, "integrate_get_samples: Insufficient memory.",
458                      TCL_STATIC);                      TCL_STATIC);
459        return TCL_ERROR;        return TCL_ERROR;
460      }      }
461    
462      /* Get the display units at string... */      /* Get the display units at string... */
463      ustring = Asc_UnitDimString(dp,0);      ustring = Asc_UnitDimString(dp,0);
464    
465      /* Get the conversion factor... */      /* Get the conversion factor... */
466      du = (struct Units *)LookupUnits(ustring);      du = (struct Units *)LookupUnits(ustring);
467      if (du == NULL) {      if (du == NULL) {
468        ERROR_REPORTER_HERE(ASC_PROG_ERR,"LookupUnits failed :-/");        ERROR_REPORTER_HERE(ASC_PROG_ERR,"LookupUnits failed :-/");
469        stat = 1;        stat = 1;
470      }else{      }else{
471        stat = 0;        stat = 0;
472        /* fill 'uvalues' with scaled values (in output units) */        /* fill 'uvalues' with scaled values (in output units) */
473        uv = uvalues;        uv = uvalues;
474        for (i=0; i < len; i++) {        for (i=0; i < len; i++) {
475          /* convert to output units */          /* convert to output units */
476          stat = Asc_UnitConvert(du,samplelist_get(&l_samplelist,i),uv,1);          stat = Asc_UnitConvert(du,samplelist_get(&l_samplelist,i),uv,1);
477          if (stat) {          if (stat) {
478            /* any problems, just stop */            /* any problems, just stop */
479            break;            break;
480          }          }
481          uv++;          uv++;
482        }        }
483      }      }
484      if (stat) {      if (stat) {
485        /* there was a problem, so free the allocated space */        /* there was a problem, so free the allocated space */
486        ascfree(uvalues);        ascfree(uvalues);
487      }      }
488    }    }
489    
490    /* give Tcl the units string */    /* give Tcl the units string */
491    Tcl_AppendResult(interp,"{",ustring,"} {",SNULL);    Tcl_AppendResult(interp,"{",ustring,"} {",SNULL);
492    
493    /* work out what precision we want */    /* work out what precision we want */
494    prec = Asc_UnitGetCPrec();    prec = Asc_UnitGetCPrec();
495    
496    len--; /* last one is a special case...? */    len--; /* last one is a special case...? */
497    if (!trydu || stat){    if (!trydu || stat){
498      /* no unit conversion, or failed unit conversion: use the raw values */      /* no unit conversion, or failed unit conversion: use the raw values */
499      for(i=0;i<len;i++){      for(i=0;i<len;i++){
500          sprintf(sval,"%.*g ",prec,samplelist_get(&l_samplelist,i));          sprintf(sval,"%.*g ",prec,samplelist_get(&l_samplelist,i));
501          Tcl_AppendResult(interp,sval,SNULL);          Tcl_AppendResult(interp,sval,SNULL);
502      }      }
503    }else{    }else{
504      for(i=0;i<len;i++){      for(i=0;i<len;i++){
505          sprintf(sval,"%.*g ",prec,uvalues[i]);          sprintf(sval,"%.*g ",prec,uvalues[i]);
506          Tcl_AppendResult(interp,sval,SNULL);          Tcl_AppendResult(interp,sval,SNULL);
507      }      }
508      ascfree(uvalues);      ascfree(uvalues);
509    }    }
510    for (i=0; i<len; i++) {    for (i=0; i<len; i++) {
511      /* print number and a blank */      /* print number and a blank */
512      sprintf(sval,"%.*g ",prec,uvalues[i]);      sprintf(sval,"%.*g ",prec,uvalues[i]);
513      Tcl_AppendResult(interp,sval,SNULL);      Tcl_AppendResult(interp,sval,SNULL);
514    }    }
515    
516    sprintf(sval,"%.*g",prec,uvalues[len]);    sprintf(sval,"%.*g",prec,uvalues[len]);
517    Tcl_AppendResult(interp,sval,"}",SNULL);    Tcl_AppendResult(interp,sval,"}",SNULL);
518    
519    if (trydu && !stat) {    if (trydu && !stat) {
520      ascfree(uvalues);      ascfree(uvalues);
521    }    }
522    return TCL_OK;    return TCL_OK;
523  }  }
524    
525  int Asc_IntegSetXSamplesCmd(ClientData cdata, Tcl_Interp *interp,  int Asc_IntegSetXSamplesCmd(ClientData cdata, Tcl_Interp *interp,
526                            int argc, CONST84 char *argv[])                            int argc, CONST84 char *argv[])
527  {  {
528    struct Units *du = NULL;    struct Units *du = NULL;
529    const dim_type *dp;    const dim_type *dp;
530    dim_type *mydp;    dim_type *mydp;
531    long i,len;    long i,len;
532    double *uvalues = NULL;    double *uvalues = NULL;
533    double *uv;    double *uv;
534    int stat;    int stat;
535    
536    UNUSED_PARAMETER(cdata);    UNUSED_PARAMETER(cdata);
537    
538    if (argc == 1) {    if (argc == 1) {
539      samplelist_assign(&l_samplelist,0L,NULL,NULL);      samplelist_assign(&l_samplelist,0L,NULL,NULL);
540      return TCL_OK;      return TCL_OK;
541    }    }
542    
543    if (argc <4) {    if (argc <4) {
544      Tcl_SetResult(interp,      Tcl_SetResult(interp,
545                    "Syntax: integrate_set_samples"                    "Syntax: integrate_set_samples"
546                    " <units> <value [value...] value>",                    " <units> <value [value...] value>",
547                    TCL_STATIC);                    TCL_STATIC);
548      FPRINTF(ASCERR,"ERROR: integrate_set_samples needs at least 3 args.");      FPRINTF(ASCERR,"ERROR: integrate_set_samples needs at least 3 args.");
549      return TCL_ERROR;      return TCL_ERROR;
550    }    }
551    
552    du = (struct Units *)LookupUnits(argv[1]);    du = (struct Units *)LookupUnits(argv[1]);
553    if (du == NULL) {    if (du == NULL) {
554      Tcl_SetResult(interp, "integrate_set_samples: first arg not valid units.",      Tcl_SetResult(interp, "integrate_set_samples: first arg not valid units.",
555                    TCL_STATIC);                    TCL_STATIC);
556      return TCL_ERROR;      return TCL_ERROR;
557    }    }
558    dp = (const dim_type *)UnitsDimensions(du);    dp = (const dim_type *)UnitsDimensions(du);
559  #if INTDEBUG  #if INTDEBUG
560    FPRINTF(ASCERR,"user dimen looks like\n");    FPRINTF(ASCERR,"user dimen looks like\n");
561    PrintDimen(ASCERR,dp);    PrintDimen(ASCERR,dp);
562  #endif  #endif
563    mydp = ASC_NEW(dim_type); /* we are giving away */    mydp = ASC_NEW(dim_type); /* we are giving away */
564    if (mydp == NULL) {    if (mydp == NULL) {
565      Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory",      Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory",
566                    TCL_STATIC);                    TCL_STATIC);
567      return TCL_ERROR;      return TCL_ERROR;
568    }    }
569    CopyDimensions(dp,mydp);    CopyDimensions(dp,mydp);
570  #if INTDEBUG  #if INTDEBUG
571    FPRINTF(ASCERR,"copy of user dimen looks like\n");    FPRINTF(ASCERR,"copy of user dimen looks like\n");
572    PrintDimen(ASCERR,mydp);    PrintDimen(ASCERR,mydp);
573  #endif  #endif
574    
575    len = argc -2;    len = argc -2;
576    uvalues = ASC_NEW_ARRAY(double,len); /* we are giving away */    uvalues = ASC_NEW_ARRAY(double,len); /* we are giving away */
577    if (uvalues==NULL) {    if (uvalues==NULL) {
578      Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory",      Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory",
579                    TCL_STATIC);                    TCL_STATIC);
580      ascfree(mydp);      ascfree(mydp);
581      return TCL_ERROR;      return TCL_ERROR;
582    }    }
583    stat = 0;    stat = 0;
584    uv = uvalues;    uv = uvalues;
585    for (i=0; i<len; i++) {    for (i=0; i<len; i++) {
586      if(Tcl_GetDouble(interp,argv[i+2],uv)!=TCL_OK) {      if(Tcl_GetDouble(interp,argv[i+2],uv)!=TCL_OK) {
587        stat = 1;        stat = 1;
588        break;        break;
589      }      }
590      stat = Asc_UnitConvert(du,*uv,uv,0);      stat = Asc_UnitConvert(du,*uv,uv,0);
591      if (stat) {      if (stat) {
592        break;        break;
593      }      }
594      uv++;      uv++;
595    }    }
596    Tcl_ResetResult(interp);    Tcl_ResetResult(interp);
597    if (stat) {    if (stat) {
598      Tcl_SetResult(interp, "integrate_set_samples: Invalid value given. (",      Tcl_SetResult(interp, "integrate_set_samples: Invalid value given. (",
599                    TCL_STATIC);                    TCL_STATIC);
600      Tcl_AppendResult(interp,argv[i+2],")",SNULL);      Tcl_AppendResult(interp,argv[i+2],")",SNULL);
601      ascfree(uvalues);      ascfree(uvalues);
602      ascfree(mydp);      ascfree(mydp);
603      return TCL_ERROR;      return TCL_ERROR;
604    }    }
605    if(!samplelist_assign(&l_samplelist,len,uvalues,mydp)){    if(!samplelist_assign(&l_samplelist,len,uvalues,mydp)){
606      Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory.",      Tcl_SetResult(interp, "integrate_set_samples: Insufficient memory.",
607                    TCL_STATIC);                    TCL_STATIC);
608      ascfree(uvalues);      ascfree(uvalues);
609      ascfree(mydp);      ascfree(mydp);
610      return TCL_ERROR;      return TCL_ERROR;
611    }    }
612    return TCL_OK;    return TCL_OK;
613  }  }
614    
615  /*----------------------------------------------------------------  /*----------------------------------------------------------------
616    FUNCTIONS THAT QUERY THE INTEGRATOR/SOLVER    FUNCTIONS THAT QUERY THE INTEGRATOR/SOLVER
617  */  */
618    
619  /**  /**
620      Tcl/Tk interface function: is the specified instance integrable.      Tcl/Tk interface function: is the specified instance integrable.
621    
622      There is a problem with this part of the interface. The      There is a problem with this part of the interface. The
623      'integrator_isintegrable' function is being called before the      'integrator_isintegrable' function is being called before the
624      'integrator_analyse' step, which means that it really doesn't have all the      'integrator_analyse' step, which means that it really doesn't have all the
625      information that it needs to be able to say. I've reworked it so that      information that it needs to be able to say. I've reworked it so that
626      you would expect to get errors back from the 'integrator_solve' step if      you would expect to get errors back from the 'integrator_solve' step if
627      there is any problem.      there is any problem.
628    
629      All this function now does is to check that the instance is OK, and to      All this function now does is to check that the instance is OK, and to
630      check that a valid integrator engine is specified.      check that a valid integrator engine is specified.
631  */  */
632  int Asc_IntegInstIntegrableCmd(ClientData cdata,Tcl_Interp *interp,  int Asc_IntegInstIntegrableCmd(ClientData cdata,Tcl_Interp *interp,
633                               int argc, CONST84 char *argv[])                               int argc, CONST84 char *argv[])
634  {  {
635    struct Instance *i=NULL;    struct Instance *i=NULL;
636    IntegratorEngine integrator=INTEG_UNKNOWN;    IntegratorEngine integrator=INTEG_UNKNOWN;
637    int result=0;         /* 0 = FALSE; 1 = TRUE */    int result=0;         /* 0 = FALSE; 1 = TRUE */
638    
639    UNUSED_PARAMETER(cdata);    UNUSED_PARAMETER(cdata);
640    
641    if ( argc != 3 ) {    if ( argc != 3 ) {
642      Tcl_SetResult(interp, "integrate_able <solver,current,search> <lsode>",      Tcl_SetResult(interp, "integrate_able <solver,current,search> <lsode>",
643                    TCL_STATIC);                    TCL_STATIC);
644      return TCL_ERROR;      return TCL_ERROR;
645    }    }
646    
647    if (strncmp(argv[1],"solver",3)==0) {    if (strncmp(argv[1],"solver",3)==0) {
648      i = g_solvinst_cur;      i = g_solvinst_cur;
649    } else {    } else {
650      if (strncmp(argv[1],"search",3)==0) {      if (strncmp(argv[1],"search",3)==0) {
651        i = g_search_inst;        i = g_search_inst;
652      } else {      } else {
653        if (strncmp(argv[1],"current",3)==0) {        if (strncmp(argv[1],"current",3)==0) {
654          i = g_curinst;          i = g_curinst;
655        } else {        } else {
656          Tcl_SetResult(interp,          Tcl_SetResult(interp,
657                        "integrate_able: arg 1 is current, search, or solver",                        "integrate_able: arg 1 is current, search, or solver",
658                        TCL_STATIC);                        TCL_STATIC);
659          return TCL_ERROR;          return TCL_ERROR;
660        }        }
661      }      }
662    }    }
663    
664    if (!i) {    if (!i) {
665      Tcl_SetResult(interp, "0", TCL_STATIC);      Tcl_SetResult(interp, "0", TCL_STATIC);
666      FPRINTF(ASCERR,"NULL instance sent to integrate_able.\n");      FPRINTF(ASCERR,"NULL instance sent to integrate_able.\n");
667      return TCL_OK;      return TCL_OK;
668    }    }
669    
670    integrator = INTEG_UNKNOWN;    integrator = INTEG_UNKNOWN;
671    
672    if (strncmp(argv[2],"blsode",3)==0) {    if (strncmp(argv[2],"blsode",3)==0) {
673      integrator = INTEG_LSODE;      integrator = INTEG_LSODE;
674    }else if (strncmp(argv[2],"ida",3)==0) {    }else if (strncmp(argv[2],"ida",3)==0) {
675      integrator = INTEG_IDA;      integrator = INTEG_IDA;
676    }    }
677    
678    result = (integrator != INTEG_UNKNOWN);    result = (integrator != INTEG_UNKNOWN);
679    if (result) {    if (result) {
680      Tcl_SetResult(interp, "1", TCL_STATIC);      Tcl_SetResult(interp, "1", TCL_STATIC);
681    } else {    } else {
682      Tcl_SetResult(interp, "0", TCL_STATIC);      Tcl_SetResult(interp, "0", TCL_STATIC);
683    }    }
684    return TCL_OK;    return TCL_OK;
685  }  }
686    
687  /*-----------------------------------*/  /*-----------------------------------*/
688    
689  /**  /**
690      Set up the Integrator.      Set up the Integrator.
691    
692      switches (in Tcl/Tk)      switches (in Tcl/Tk)
693          -engine $name          -engine $name
694          -i0 $stepindex          -i0 $stepindex
695          -i1 $stepindex          -i1 $stepindex
696          -dt0 $initstepsize          -dt0 $initstepsize
697          -dtmin $minstep          -dtmin $minstep
698          -dtmax $maxstep          -dtmax $maxstep
699          -moststeps $moststeps          -moststeps $moststeps
700  */  */
701  int Asc_IntegSetupCmd(ClientData cdata,Tcl_Interp *interp,  int Asc_IntegSetupCmd(ClientData cdata,Tcl_Interp *interp,
702                     int argc, CONST84 char *argv[])                     int argc, CONST84 char *argv[])
703  {  {
704    IntegratorEngine integrator = INTEG_UNKNOWN;    IntegratorEngine integrator = INTEG_UNKNOWN;
705    char buf[MAXIMUM_NUMERIC_LENGTH];         /* string to hold integer */    char buf[MAXIMUM_NUMERIC_LENGTH];         /* string to hold integer */
706    CONST84 char *engine = NULL;    CONST84 char *engine = NULL;
707    int result = 0;         /* 0 = FALSE; 1 = TRUE */    int result = 0;         /* 0 = FALSE; 1 = TRUE */
708    long i0=(-1), i1=(-1);    long i0=(-1), i1=(-1);
709    int ifound = 0;    int ifound = 0;
710    int k;    int k;
711    int moststeps=0;    int moststeps=0;
712    double dt0=0, dtmin=0, dtmax=0;    double dt0=0, dtmin=0, dtmax=0;
713    CONST84 char *cdt0=NULL, *cdtmin=NULL, *cdtmax=NULL, *cmoststeps=NULL,    CONST84 char *cdt0=NULL, *cdtmin=NULL, *cdtmax=NULL, *cmoststeps=NULL,
714         *ci0=NULL, *ci1=NULL;         *ci0=NULL, *ci1=NULL;
715    IntegratorReporter *reporter;    IntegratorReporter *reporter;
716    IntegratorSystem *blsys;    IntegratorSystem *blsys;
717    
718    UNUSED_PARAMETER(cdata);    UNUSED_PARAMETER(cdata);
719    
720    k = 1;    k = 1;
721    while (k < (argc-1)) { /* arguments come in pairs */    while (k < (argc-1)) { /* arguments come in pairs */
722      if (strcmp(argv[k],"-engine")==0) {      if (strcmp(argv[k],"-engine")==0) {
723        engine = argv[k+1];        engine = argv[k+1];
724        k+=2;        k+=2;
725        continue;        continue;
726      }      }
727      if (strcmp(argv[k],"-i1")==0) {      if (strcmp(argv[k],"-i1")==0) {
728        ci1 = argv[k+1];        ci1 = argv[k+1];
729        k+=2;        k+=2;
730        continue;        continue;
731      }      }
732      if (strcmp(argv[k],"-i0")==0) {      if (strcmp(argv[k],"-i0")==0) {
733        ci0 = argv[k+1];        ci0 = argv[k+1];
734        k+=2;        k+=2;
735        continue;        continue;
736      }      }
737      if (strcmp(argv[k],"-moststeps")==0) {      if (strcmp(argv[k],"-moststeps")==0) {
738        cmoststeps = argv[k+1];        cmoststeps = argv[k+1];
739        k+=2;        k+=2;
740        continue;        continue;
741      }      }
742      if (strcmp(argv[k],"-dtmax")==0) {      if (strcmp(argv[k],"-dtmax")==0) {
743        cdtmax = argv[k+1];        cdtmax = argv[k+1];
744        k+=2;        k+=2;
745        continue;        continue;
746      }      }
747      if (strcmp(argv[k],"-dtmin")==0) {      if (strcmp(argv[k],"-dtmin")==0) {
748        cdtmin = argv[k+1];        cdtmin = argv[k+1];
749        k+=2;        k+=2;
750        continue;        continue;
751      }      }
752      if (strcmp(argv[k],"-dt0")==0) {      if (strcmp(argv[k],"-dt0")==0) {
753        cdt0 = argv[k+1];        cdt0 = argv[k+1];
754        k+=2;        k+=2;
755        continue;        continue;
756      }      }
757      Tcl_AppendResult(interp,argv[0],": unrecognized option ",      Tcl_AppendResult(interp,argv[0],": unrecognized option ",
758                       argv[k],".",SNULL);                       argv[k],".",SNULL);
759      return TCL_ERROR;      return TCL_ERROR;
760    }    }
761    
762    if (engine != NULL && strncmp(engine,"BLSODE",3)==0) {    if (engine != NULL && strncmp(engine,"BLSODE",3)==0) {
763      integrator = INTEG_LSODE;      integrator = INTEG_LSODE;
764      ifound=1;      ifound=1;
765    }    }
766    
767    if (engine != NULL && strncmp(engine,"IDA",3)==0) {    if (engine != NULL && strncmp(engine,"IDA",3)==0) {
768      integrator = INTEG_IDA;      integrator = INTEG_IDA;
769      ifound=1;      ifound=1;
770    }    }
771    
772    if (!ifound) {    if (!ifound) {
773      Tcl_SetResult(interp, "Unsupported integrator", TCL_STATIC);      Tcl_SetResult(interp, "Unsupported integrator", TCL_STATIC);
774      Tcl_AppendResult(interp," ",engine,SNULL);      Tcl_AppendResult(interp," ",engine,SNULL);
775      return TCL_ERROR;      return TCL_ERROR;
776    }    }
777    if (ci0 != NULL && ci1 != NULL) {    if (ci0 != NULL && ci1 != NULL) {
778      /* get i0, i1 if both supplied. */      /* get i0, i1 if both supplied. */
779      long i;      long i;
780      if (Tcl_ExprLong(interp,ci0,&i)==TCL_ERROR|| i<0) {      if (Tcl_ExprLong(interp,ci0,&i)==TCL_ERROR|| i<0) {
781        Tcl_ResetResult(interp);        Tcl_ResetResult(interp);
782        Tcl_SetResult(interp, "integrate_setup: index i0 invalid", TCL_STATIC);        Tcl_SetResult(interp, "integrate_setup: index i0 invalid", TCL_STATIC);
783        return TCL_ERROR;        return TCL_ERROR;
784      }      }
785      i0=i;      i0=i;
786      if (Tcl_ExprLong(interp,ci1,&i)==TCL_ERROR|| i<i0) {      if (Tcl_ExprLong(interp,ci1,&i)==TCL_ERROR|| i<i0) {
787        Tcl_ResetResult(interp);        Tcl_ResetResult(interp);
788        Tcl_SetResult(interp, "integrate_setup: index i1 invalid", TCL_STATIC);        Tcl_SetResult(interp, "integrate_setup: index i1 invalid", TCL_STATIC);
789        return TCL_ERROR;        return TCL_ERROR;
790      }      }
791      i1=i;      i1=i;
792    }    }
793    if (cdt0 != NULL) {    if (cdt0 != NULL) {
794      if (Tcl_GetDouble(interp,cdt0,&dt0) != TCL_OK) {      if (Tcl_GetDouble(interp,cdt0,&dt0) != TCL_OK) {
795        Tcl_ResetResult(interp);        Tcl_ResetResult(interp);
796        Tcl_AppendResult(interp, "integrate_setup: initial step length invalid",        Tcl_AppendResult(interp, "integrate_setup: initial step length invalid",
797                         " (",cdt0,")", SNULL);                         " (",cdt0,")", SNULL);
798        return TCL_ERROR;        return TCL_ERROR;
799      }      }
800    }    }
801    if (cdtmin != NULL) {    if (cdtmin != NULL) {
802      if (Tcl_GetDouble(interp,cdtmin,&dtmin) != TCL_OK || dtmin < 0) {      if (Tcl_GetDouble(interp,cdtmin,&dtmin) != TCL_OK || dtmin < 0) {
803        Tcl_ResetResult(interp);        Tcl_ResetResult(interp);
804        Tcl_AppendResult(interp, "integrate_setup: minimum step length invalid",        Tcl_AppendResult(interp, "integrate_setup: minimum step length invalid",
805                         " (",cdtmin,")", SNULL);                         " (",cdtmin,")", SNULL);
806        return TCL_ERROR;        return TCL_ERROR;
807      }      }
808    }    }
809    if (cdtmax != NULL) {    if (cdtmax != NULL) {
810      if (Tcl_GetDouble(interp,cdtmax,&dtmax) != TCL_OK || dtmax < dtmin) {      if (Tcl_GetDouble(interp,cdtmax,&dtmax) != TCL_OK || dtmax < dtmin) {
811        Tcl_ResetResult(interp);        Tcl_ResetResult(interp);
812        Tcl_AppendResult(interp, "integrate_setup: maximum step length invalid",        Tcl_AppendResult(interp, "integrate_setup: maximum step length invalid",
813                         " (",cdtmax,")", SNULL);                         " (",cdtmax,")", SNULL);
814        return TCL_ERROR;        return TCL_ERROR;
815      }      }
816    }    }
817    if (cmoststeps != NULL) {    if (cmoststeps != NULL) {
818      if (Tcl_GetInt(interp,cmoststeps,&moststeps) != TCL_OK || moststeps < 0) {      if (Tcl_GetInt(interp,cmoststeps,&moststeps) != TCL_OK || moststeps < 0) {
819        Tcl_ResetResult(interp);        Tcl_ResetResult(interp);
820        Tcl_AppendResult(interp, "integrate_setup: maximum internal steps bad",        Tcl_AppendResult(interp, "integrate_setup: maximum internal steps bad",
821                         " (",cmoststeps,")", SNULL);                         " (",cmoststeps,")", SNULL);
822        return TCL_ERROR;        return TCL_ERROR;
823      }      }
824    }    }
825    
826    reporter = Asc_GetIntegReporter();    reporter = Asc_GetIntegReporter();
827    
828    blsys = integrator_new(g_solvsys_cur,g_solvinst_cur);    blsys = integrator_new(g_solvsys_cur,g_solvinst_cur);
829    integrator_set_engine(blsys, integrator);    integrator_set_engine(blsys, integrator);
830    integrator_set_reporter(blsys, reporter);    integrator_set_reporter(blsys, reporter);
831    integrator_set_samples(blsys,&l_samplelist);    integrator_set_samples(blsys,&l_samplelist);
832    integrator_set_stepzero(blsys,dt0);    integrator_set_stepzero(blsys,dt0);
833    integrator_set_minstep(blsys,dtmin);    integrator_set_minstep(blsys,dtmin);
834    integrator_set_maxstep(blsys,dtmax);    integrator_set_maxstep(blsys,dtmax);
835    integrator_set_maxsubsteps(blsys,moststeps);    integrator_set_maxsubsteps(blsys,moststeps);
836    
837    result = integrator_analyse(blsys);    result = integrator_analyse(blsys);
838    
839    /* go and solve it */    /* go and solve it */
840    integrator_solve(blsys, i0, i1);    integrator_solve(blsys, i0, i1);
841    
842    /* once solution is finished, free whatever we allocated */    /* once solution is finished, free whatever we allocated */
843    integrator_free(blsys);    integrator_free(blsys);
844    
845    sprintf(buf, "%d", result);    sprintf(buf, "%d", result);
846    Tcl_SetResult(interp, buf, TCL_VOLATILE);    Tcl_SetResult(interp, buf, TCL_VOLATILE);
847    return TCL_OK;    return TCL_OK;
848  }  }
849    
850  /********************************************************************/  /********************************************************************/
851    
852  int Asc_IntegCleanupCmd(ClientData cdata,Tcl_Interp *interp,  int Asc_IntegCleanupCmd(ClientData cdata,Tcl_Interp *interp,
853                       int argc, CONST84 char *argv[])                       int argc, CONST84 char *argv[])
854  {  {
855    UNUSED_PARAMETER(cdata);    UNUSED_PARAMETER(cdata);
856    (void)argv;     /* stop gcc whine about unused parameter */    (void)argv;     /* stop gcc whine about unused parameter */
857    
858    if (argc!=1) {    if (argc!=1) {
859      Tcl_SetResult(interp, "integrate_cleanup takes no arguments", TCL_STATIC);      Tcl_SetResult(interp, "integrate_cleanup takes no arguments", TCL_STATIC);
860      return TCL_ERROR;      return TCL_ERROR;
861    }    }
862    
863    /* integrator_cleanup(); */    /* integrator_cleanup(); */
864    return TCL_OK;    return TCL_OK;
865  }  }
866    
867  /*-------------------------------------------------------------------  /*-------------------------------------------------------------------
868    REPORTER FUNCTIONS    REPORTER FUNCTIONS
869  */  */
870    
871  FILE *integ_y_out;  FILE *integ_y_out;
872  FILE *integ_obs_out;  FILE *integ_obs_out;
873    
874  IntegratorReporter *Asc_GetIntegReporter(){  IntegratorReporter *Asc_GetIntegReporter(){
875      IntegratorReporter *r;      IntegratorReporter *r;
876      r = (IntegratorReporter *)ascmalloc(sizeof(IntegratorReporter));      r = (IntegratorReporter *)ascmalloc(sizeof(IntegratorReporter));
877      r->init = &Asc_IntegReporterInit;      r->init = &Asc_IntegReporterInit;
878      r->write = &Asc_IntegReporterWrite;      r->write = &Asc_IntegReporterWrite;
879      r->write_obs = &Asc_IntegReporterWriteObs;      r->write_obs = &Asc_IntegReporterWriteObs;
880      r->close = &Asc_IntegReporterClose;      r->close = &Asc_IntegReporterClose;
881      CONSOLE_DEBUG("CREATED INTEGRATORREPORTER FOR TCL/TK INTERFACE");      CONSOLE_DEBUG("CREATED INTEGRATORREPORTER FOR TCL/TK INTERFACE");
882      return r;      return r;
883  }  }
884    
885  int Asc_IntegReporterInit(IntegratorSystem *blsys){  int Asc_IntegReporterInit(IntegratorSystem *blsys){
886      CONSOLE_DEBUG("INITIALISING REPORTER");      CONSOLE_DEBUG("INITIALISING REPORTER");
887    
888      /* set up output files */      /* set up output files */
889      integ_y_out = Asc_IntegOpenYFile();      integ_y_out = Asc_IntegOpenYFile();
890      integ_obs_out = Asc_IntegOpenObsFile();      integ_obs_out = Asc_IntegOpenObsFile();
891    
892      CONSOLE_DEBUG("RELEASING FILES");      CONSOLE_DEBUG("RELEASING FILES");
893    
894      Asc_IntegReleaseYFile();      Asc_IntegReleaseYFile();
895      Asc_IntegReleaseObsFile();      Asc_IntegReleaseObsFile();
896    
897      CONSOLE_DEBUG("WRITING HEADERS");      CONSOLE_DEBUG("WRITING HEADERS");
898    
899      /* write headers to yout, obsout and initial points */      /* write headers to yout, obsout and initial points */
900      Asc_IntegPrintYHeader(integ_y_out,blsys);      Asc_IntegPrintYHeader(integ_y_out,blsys);
901      Asc_IntegPrintYLine(integ_y_out,blsys);      Asc_IntegPrintYLine(integ_y_out,blsys);
902      Asc_IntegPrintObsHeader(integ_obs_out,blsys);      Asc_IntegPrintObsHeader(integ_obs_out,blsys);
903      Asc_IntegPrintObsLine(integ_obs_out,blsys);      Asc_IntegPrintObsLine(integ_obs_out,blsys);
904  }  }
905    
906  int Asc_IntegReporterWrite(IntegratorSystem *blsys){  int Asc_IntegReporterWrite(IntegratorSystem *blsys){
907      /* write out a line of stuff */      /* write out a line of stuff */
908      Asc_IntegPrintYLine(integ_y_out,blsys);      Asc_IntegPrintYLine(integ_y_out,blsys);
909  }  }
910    
911  int Asc_IntegReporterWriteObs(IntegratorSystem *blsys){  int Asc_IntegReporterWriteObs(IntegratorSystem *blsys){
912      Asc_IntegPrintObsLine(integ_obs_out,blsys);      Asc_IntegPrintObsLine(integ_obs_out,blsys);
913  }  }
914    
915  int Asc_IntegReporterClose(IntegratorSystem *blsys){  int Asc_IntegReporterClose(IntegratorSystem *blsys){
916      /* close the file streams */      /* close the file streams */
917      if (integ_y_out!=NULL) {      if (integ_y_out!=NULL) {
918          fclose(integ_y_out);          fclose(integ_y_out);
919      }      }
920    
921      if (integ_obs_out!=NULL) {      if (integ_obs_out!=NULL) {
922          fclose(integ_obs_out);          fclose(integ_obs_out);
923      }      }
924  }  }
925    

Legend:
Removed from v.689  
changed lines
  Added in v.690

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22