/[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 690 - (show annotations) (download) (as text)
Thu Jun 22 00:48:31 2006 UTC (14 years, 5 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 /* 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