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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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