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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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