root/branches/newmole/source/assert_results.cpp

Revision 2542, 85.2 kB (checked in by rjrw, 3 weeks ago)

Merged from trunk r2448:2541

  • Property svn:eol-style set to native
Line 
1/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2 * others.  For conditions of distribution and use see copyright notice in license.txt */
3/*InitAssertResults, this must be called first, done at startup of ParseCommands*/
4/*ParseAssertResults - parse input stream */
5/*lgCheckAsserts, checks asserts, last thing cloudy calls, returns true if all are ok, false if problems */
6#include "cddefines.h"
7#include "input.h"
8#include "conv.h"
9#include "optimize.h"
10#include "iso.h"
11#include "called.h"
12#include "atmdat.h"
13#include "hcmap.h"
14#include "thermal.h"
15#include "pressure.h"
16#include "struc.h"
17#include "wind.h"
18#include "h2.h"
19#include "colden.h"
20#include "dense.h"
21#include "lines.h"
22#include "secondaries.h"
23#include "radius.h"
24#include "version.h"
25#include "hmi.h"
26#include "prt.h"
27#include "grainvar.h"
28#include "atomfeii.h"
29#include "cddrive.h"
30#include "elementnames.h"
31#include "assertresults.h"
32#include "taulines.h"
33#include "lines_service.h"
34
35/* flag to remember that InitAssertResults was called */
36static bool lgInitDone=false ,
37        /* will be set true when space for asserts is allocated */
38        lgSpaceAllocated=false;
39
40/* number of asserts we can handle, used in dim of space */
41#define NASSERTS 100
42
43/* default relative error for asserted quantities */
44static realnum ErrorDefault;
45
46/* dim of 5 also appears in MALLOC below */
47#define NCHAR 5
48static char **chAssertLineLabel;
49
50/* this will be = for equal, < or > for limit */
51static char *chAssertLimit;
52
53/* this will be a two character label identifying which type of assert */
54static char **chAssertType;
55
56/* the values and error in the asserted quantity */
57static double *AssertQuantity ,*AssertQuantity2 ,*AssertError, **Param;
58
59/* this flag says where we print linear or log quantity */
60static int *lgQuantityLog;
61static long nAsserts=0;
62static realnum *wavelength;
63
64/*======================================================================*/
65/*InitAssertResults, this must be called first, done at startup of ParseCommands*/
66void InitAssertResults(void)
67{
68        /* set flag that init was called, and set number of asserts to zero.
69         * this is done by ParseComments for every model, even when no asserts will
70         * be done, so do not allocate space at this time */
71        lgInitDone = true;
72        nAsserts = 0;
73
74        /* default error, changed with ASSERT SET ERROR */
75        ErrorDefault = 0.05f;
76}
77
78/*======================================================================*/
79/*ParseAssertResults parse the assert command */
80void ParseAssertResults(void)
81{
82        long i ,
83                nelem,
84                n2;
85        bool lgEOL;
86        char chLabel[INPUT_LINE_LENGTH];
87
88        DEBUG_ENTRY( "ParseAssertResults()" );
89
90        if( !lgInitDone )
91        {
92                fprintf( ioQQQ, " ParseAssertResults called before InitAsserResults\n" );
93                cdEXIT(EXIT_FAILURE);
94        }
95
96        /* has space been allocated yet? */
97        if( !lgSpaceAllocated )
98        {
99                /* - no, we must allocate it */
100                /* remember that space has been allocated */
101                lgSpaceAllocated = true;
102
103                /* create space for the array of labels*/
104                chAssertLineLabel = ((char **)MALLOC(NASSERTS*sizeof(char *)));
105
106                /* the 2-character string saying what type of assert */
107                chAssertType = ((char **)MALLOC(NASSERTS*sizeof(char *)));
108
109                /* these are a pair of optional parameters */
110                Param = ((double **)MALLOC(NASSERTS*sizeof(double * )));
111
112                /* now fill out the 2D arrays */
113                for( i=0; i<NASSERTS; ++i )
114                {
115                        chAssertLineLabel[i] = ((char *)MALLOC(NCHAR*sizeof(char )));
116                        strcpy(chAssertLineLabel[i],"unkn" );
117
118                        /* two char plus eol */
119                        chAssertType[i] = ((char *)MALLOC(3*sizeof(char )));
120                        strcpy(chAssertType[i],"un" );
121
122                        /* these are a pair of optional parameters */
123                        Param[i] = ((double *)MALLOC(5*sizeof(double )));
124                }
125
126                /* now make space for the asserted quantities  */
127                AssertQuantity = ((double *)MALLOC(NASSERTS*sizeof(double )));
128
129                AssertQuantity2 = ((double *)MALLOC(NASSERTS*sizeof(double )));
130
131                /* now the errors */
132                AssertError = ((double *)MALLOC(NASSERTS*sizeof(double )));
133
134                /* now the line wavelengths */
135                wavelength = ((realnum *)MALLOC(NASSERTS*sizeof(realnum )));
136
137                /* now the flag saying whether should be log */
138                lgQuantityLog = ((int *)MALLOC(NASSERTS*sizeof(int )));
139
140                /* the flag for upper, lower limit, or equal */
141                chAssertLimit = ((char *)MALLOC(NASSERTS*sizeof(char )));
142        }
143        /* end space allocation - we are ready to roll */
144
145        /* first say we do not know what job to do */
146        strcpy( chAssertType[nAsserts] , "un" );
147
148        /* false means print linear quantity - will be set true if entered
149         * quantity comes in as a log */
150        lgQuantityLog[nAsserts] = false;
151
152        /* all asserts have option for quantity to be a limit, or the quantity itself */
153        if( nMatch("<",input.chCARDCAPS ) )
154        {
155                chAssertLimit[nAsserts] = '<';
156        }
157        else if( nMatch(">",input.chCARDCAPS ) )
158        {
159                chAssertLimit[nAsserts] = '>';
160        }
161        else
162        {
163                chAssertLimit[nAsserts] = '=';
164        }
165
166        /* which quantity will we check?, first is */
167
168        if( nMatch(" SET",input.chCARDCAPS ) )
169        {
170                /* set an option for the assert command, not an actual asserted
171                 * quantity */
172
173                /* decrement number of asserts since will be incremented below,
174                 * this is not an actual asserted quantity */
175                if( nAsserts >0 )
176                        fprintf(ioQQQ," The default assert error is being changed after"
177                        " some asserts were entered.  \n This will only affect asserts"
178                        " that come after this command.\n");
179                --nAsserts;
180
181                if( nMatch("ERRO",input.chCARDCAPS ) )
182                {
183                        /* set default error */
184                        i = 5;
185                        ErrorDefault =
186                                (realnum)FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
187                        if( lgEOL )
188                                NoNumb(input.chOrgCard);
189                }
190                else
191                {
192                        /* problem - no recognized quantity */
193                        fprintf( ioQQQ,
194                                " I could not identify an option on this ASSERT SET XXX command.\n");
195                        fprintf( ioQQQ, " Sorry.\n" );
196                        cdEXIT(EXIT_FAILURE);
197                }
198        }
199
200        /* assert mean ionization */
201        else if( nMatch("IONI",input.chCARDCAPS ) )
202        {
203
204                /* say that this will be mean ionization fraction */
205
206                /* f will indicate average over radius, F over volume -
207                 * check whether keyword radius or volume occurs,
208                 * default will be radius */
209                if( nMatch("VOLU",input.chCARDCAPS ) )
210                {
211                        strcpy( chAssertType[nAsserts] , "F " );
212                }
213                else
214                {
215                        /* this is default case, Fraction over radius */
216                        strcpy( chAssertType[nAsserts] , "f " );
217                }
218
219                /* first get element label and make null terminated string*/
220                if( (nelem = GetElem(input.chCARDCAPS)) < 0 )
221                {
222                        fprintf( ioQQQ,
223                                " I could not identify an element name on this line.\n");
224                        fprintf( ioQQQ, " Sorry.\n" );
225                        cdEXIT(EXIT_FAILURE);
226                }
227                ASSERT( nelem>= 0);
228                ASSERT( nAsserts>= 0);
229                /* we now have element name, copy 4-char string (elementnames.chElementNameShort[nelem])
230                 * into array that will be used to get ionization after calculation */
231                strcpy( chAssertLineLabel[nAsserts], elementnames.chElementNameShort[nelem] );
232
233                /* now get ionization stage, which will be saved into wavelength */
234                i = 5;
235                wavelength[nAsserts] =
236                        (realnum)FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
237                if( lgEOL )
238                {
239                        NoNumb(input.chOrgCard);
240                }
241                /* ionization stage must be 1 or greater, but not greater than nelem (c scale)+2 */
242                if( wavelength[nAsserts] < 1 || wavelength[nAsserts] > nelem+2 )
243                {
244                        fprintf( ioQQQ,
245                                "  The ionization stage is inappropriate for this element.\n");
246                        fprintf( ioQQQ, " Sorry.\n" );
247                        cdEXIT(EXIT_FAILURE);
248                }
249
250                /* now get ionization fraction, log if number is negative or == 0,
251                 * linear if positive but less than or equal to 1.*/
252                AssertQuantity[nAsserts] =
253                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
254                if( lgEOL )
255                        NoNumb(input.chOrgCard);
256
257                /* optional error, default available (cannot do before loop since we
258                * do not know how many numbers are on line */
259                AssertError[nAsserts] =
260                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
261                if( lgEOL )
262                        /* default error was set in define above */
263                        AssertError[nAsserts] = ErrorDefault;
264
265                if( nMatch( "GRID" , input.chCARDCAPS ) )
266                {
267                        long int j;
268                        /* skip nOptimiz values */
269                        for( j=0; j<optimize.nOptimiz; ++j )
270                        {
271                                AssertQuantity[nAsserts] =
272                                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
273                        }
274                        /*fprintf(ioQQQ,"DEBUG grid HIT %li val %.2f\n",
275                        optimize.nOptimiz , AssertQuantity[nAsserts] );*/
276                        if( lgEOL )
277                        {
278                                fprintf(ioQQQ,"PROBLEM this assert grid command does not have enough values.  sorry\n");
279                        }
280                }
281
282                /* now make sure we end up with positive linear ionization fraction that
283                 * is greater than 0 but less than or equal to 1. */
284                if( AssertQuantity[nAsserts] <= 0. )
285                {
286                        /* log since negative or 0 */
287                        AssertQuantity[nAsserts] =
288                                pow(10.,AssertQuantity[nAsserts] );
289                        /* entered as a log, so print as a log too */
290                        lgQuantityLog[nAsserts] = true;
291                }
292                else if( AssertQuantity[nAsserts] > 1. )
293                {
294                        fprintf( ioQQQ,
295                                "  The ionization fraction must be less than one.\n");
296                        fprintf( ioQQQ, " Sorry.\n" );
297                        cdEXIT(EXIT_FAILURE);
298                }
299
300                /* result cannot be zero */
301                if( fabs(AssertQuantity[nAsserts]) <= SMALLDOUBLE )
302                {
303                        fprintf( ioQQQ,
304                                "  The ionization ionization fraction is too small, or zero.  Check input\n");
305                        fprintf( ioQQQ, " Sorry.\n" );
306                        cdEXIT(EXIT_FAILURE);
307                }
308        }
309
310        /* molecular fraction averaged over model */
311        else if( nMatch("MOLE",input.chCARDCAPS )&&nMatch("FRAC",input.chCARDCAPS ) )
312        {
313
314                /* say that this will be mean molecular fraction */
315
316                /* mf will indicate average over radius, MF over vol -
317                 * check whether keyword radius or volume occurs,
318                 * default will be radius */
319                /** \todo       2       NB this is not used, should do both, and more molecules (H2 only for now) */
320                if( nMatch("VOLU",input.chCARDCAPS ) )
321                {
322                        strcpy( chAssertType[nAsserts] , "MF" );
323                }
324                else
325                {
326                        /* this is default case, Fraction over radius */
327                        strcpy( chAssertType[nAsserts] , "mf" );
328                }
329
330                i = nMatch(" H2 ",input.chCARDCAPS );
331                if( i )
332                {
333                        strcpy( chAssertLineLabel[nAsserts], "H2  " );
334                        /* increment to get past the label */
335                        i += 3;
336                }
337                else if( nMatch(" CO ",input.chCARDCAPS ) )
338                {
339                        strcpy( chAssertLineLabel[nAsserts], "CO  " );
340                        i = 5;
341                }
342                else
343                {
344                        fprintf( ioQQQ,
345                                " I could not identify CO or H2 on this line.\n");
346                        fprintf( ioQQQ, " Sorry.\n" );
347                        cdEXIT(EXIT_FAILURE);
348                }
349
350                /* not meaningful */
351                wavelength[nAsserts] = 0;
352
353                /* i was set above for start of scan */
354                /* now get log of molecular fraction */
355                AssertQuantity[nAsserts] =
356                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
357                if( lgEOL )
358                {
359                        NoNumb(input.chOrgCard);
360                }
361                if( AssertQuantity[nAsserts] <= 0. )
362                {
363                        /* if negative then entered as log, but we will compare with linear */
364                        AssertQuantity[nAsserts] =
365                                pow(10.,AssertQuantity[nAsserts] );
366                }
367
368                /* optional error, default available (cannot do before loop since we
369                 * do not know how many numbers are on line */
370                AssertError[nAsserts] =
371                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
372                if( lgEOL )
373                {
374                        /* default error was set in define above */
375                        AssertError[nAsserts] = ErrorDefault;
376                }
377                /* print results as logs */
378                lgQuantityLog[nAsserts] = true;
379        }
380
381        /* assert line "LINE"  --  key is ine_ since linear option appears on some commands */
382        else if( nMatch(" LINE ",input.chCARDCAPS ) )
383        {
384                if(  nMatch("LUMI",input.chCARDCAPS ) || nMatch("INTE",input.chCARDCAPS ))
385                {
386                        /* say that this is a line luminosity or intensity*/
387                        strcpy( chAssertType[nAsserts] , "Ll" );
388
389                        /* entered as a log, so print as a log too */
390                        lgQuantityLog[nAsserts] = true;
391                }
392                else
393                {
394                        /* say that this is line relative to norm line - this is the default */
395                        strcpy( chAssertType[nAsserts] , "Lr" );
396
397                        /* entered linear quantity, so print as linear too */
398                        lgQuantityLog[nAsserts] = false;
399                }
400
401                /* this will check a line intensity, first get labels within quotes,
402                 * will be null terminated */
403                GetQuote( chLabel , input.chCARDCAPS , true );
404
405                /* check that there were exactly 4 characters in the label*/
406                if( chLabel[4] != '\0' )
407                {
408                        fprintf( ioQQQ,
409                                " The label must be exactly 4 char long, between double quotes.\n");
410                        fprintf( ioQQQ, " Sorry.\n" );
411                        cdEXIT(EXIT_FAILURE);
412                }
413
414                /* copy string into array */
415                strcpy( chAssertLineLabel[nAsserts], chLabel );
416
417                /* now get line wavelength */
418                i = 5;
419                wavelength[nAsserts] =
420                        (realnum)FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
421                if( lgEOL )
422                {
423                        NoNumb(input.chOrgCard);
424                }
425                /* check for optional micron or cm units, else interpret as Angstroms */
426                if( input.chCARDCAPS[i-1] == 'M' )
427                {
428                        /* microns */
429                        wavelength[nAsserts] *= 1e4;
430                }
431                else if( input.chCARDCAPS[i-1] == 'C' )
432                {
433                        /* centimeters */
434                        wavelength[nAsserts] *= 1e8;
435                }
436
437                /* now get intensity or luminosity -
438                 * rel intensity is linear and intensity or luminosity are log */
439                AssertQuantity[nAsserts] =
440                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
441                if( lgEOL )
442                {
443                        NoNumb(input.chOrgCard);
444                }
445                /* luminosity was entered as a log */
446                if( lgQuantityLog[nAsserts] )
447                {
448                        if( AssertQuantity[nAsserts] > DBL_MAX_10_EXP ||
449                                AssertQuantity[nAsserts] < -DBL_MAX_10_EXP )
450                        {
451                                fprintf( ioQQQ,
452                                        " The asserted quantity is a log, but is too large or small, value is %e.\n",AssertQuantity[nAsserts] );
453                                fprintf( ioQQQ, " I would crash if I tried to use it.\n Sorry.\n" );
454                                cdEXIT(EXIT_FAILURE);
455                        }
456                        AssertQuantity[nAsserts] =
457                                pow(10.,AssertQuantity[nAsserts] );
458                }
459                if( AssertQuantity[nAsserts]<= 0. )
460                {
461                        fprintf( ioQQQ,
462                                " The relative intensity must be positive, and was %e.\n",AssertQuantity[nAsserts] );
463                        fprintf( ioQQQ, " Sorry.\n" );
464                        cdEXIT(EXIT_FAILURE);
465                }
466
467                /* optional error, default available */
468                AssertError[nAsserts] =
469                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
470                if( lgEOL )
471                {
472                        /* default error was set in define above */
473                        AssertError[nAsserts] = ErrorDefault;
474                }
475        }
476
477        /* assert departure coefficients */
478        else if( nMatch("CASE",input.chCARDCAPS ) )
479        {
480                /* this is Case B for some element */
481                strcpy( chAssertType[nAsserts] , "CB" );
482                i = 5;
483                /* this is relative error */
484                AssertError[nAsserts] =
485                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
486                if( lgEOL )
487                        /* default error was set in define above */
488                        AssertError[nAsserts] = ErrorDefault;
489                AssertQuantity[nAsserts] = 0;
490                wavelength[nAsserts] = 0.;
491
492                /* faint option - do not test line if relative intensity is less
493                 * than entered value */
494                if( (i=nMatch("FAINT",input.chCARDCAPS ))>0 )
495                {
496                        Param[nAsserts][4] =
497                                FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
498                        if( lgEOL )
499                        {
500                                /* did not get 2 numbers */
501                                fprintf(ioQQQ," The assert Case B faint option must have a number,"
502                                        " the relative intensity of the fainest line to assert.\n");
503                                cdEXIT(EXIT_FAILURE);
504                        }
505                        /* number is log if <= 0 */
506                        if( Param[nAsserts][4]<=0. )
507                                Param[nAsserts][4] = pow(10., Param[nAsserts][4] );
508                }
509                else
510                {
511                        /* use default - include everything*/
512                        Param[nAsserts][4] = SMALLFLOAT;
513                }
514
515                /* range option - to limit check on a certain wavelength range */
516                if( (i=nMatch("RANG",input.chCARDCAPS ))>0 )
517                {
518                        Param[nAsserts][2] =
519                                FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
520                        Param[nAsserts][3] =
521                                FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
522                        if( lgEOL )
523                        {
524                                /* did not get 2 numbers */
525                                fprintf(ioQQQ," The assert Case B range option must have two numbers,"
526                                        " the lower and upper limit to the wavelengths in Angstroms.\n");
527                                fprintf(ioQQQ," There must be a total of three numbers on the line,"
528                                        " the relative error followed by the lower and upper limits to the "
529                                        "wavelength in Angstroms.\n");
530                                cdEXIT(EXIT_FAILURE);
531                        }
532                        if( Param[nAsserts][2]>Param[nAsserts][3])
533                        {
534                                /* make sure in increasing order */
535                                double sav = Param[nAsserts][3];
536                                Param[nAsserts][3] = Param[nAsserts][2];
537                                Param[nAsserts][2] = sav;
538                        }
539                }
540                else
541                {
542                        /* use default - include everything*/
543                        Param[nAsserts][2] = 0.;
544                        Param[nAsserts][3] = 1e30;
545                }
546                /* assert case b for H - O checking against Hummer & Storey tables */
547                if( nMatch("H-LI",input.chCARDCAPS ) )
548                {
549                        /* H-like - now get an element */
550                        if( (nelem = GetElem( input.chCARDCAPS )) < 0 )
551                        {
552                                /* no name found */
553                                fprintf(ioQQQ, "assert H-like case B did not find an element on this line, sorry %s\n",
554                                        input.chCARDCAPS );
555                                cdEXIT(EXIT_FAILURE);
556                        }
557                        if( nelem>7 )
558                        {
559                                /* beyond reach of tables */
560                                fprintf(ioQQQ, "assert H-like cannot do elements heavier than O, sorry %s\n",
561                                        input.chCARDCAPS );
562                                cdEXIT(EXIT_FAILURE);
563                        }
564                        Param[nAsserts][0] = ipH_LIKE;
565                        Param[nAsserts][1] = nelem;
566                        /* generate string to find simple prediction, as in "Ca B" */
567                        strcpy( chAssertLineLabel[nAsserts], "Ca ");
568                        if( nMatch("CASE A ",input.chCARDCAPS ) )
569                                strcat( chAssertLineLabel[nAsserts] , "A");
570                        else
571                                strcat( chAssertLineLabel[nAsserts] , "B");
572                }
573                else if( nMatch("HE-L",input.chCARDCAPS ) )
574                {
575                        /* He-like - only helium itself */
576                        Param[nAsserts][0] = ipHE_LIKE;
577                        Param[nAsserts][1] = ipHELIUM;
578
579                        strcpy( chAssertLineLabel[nAsserts] , "+Col");
580                }
581                else
582                {
583                        /*no option found */
584                        fprintf( ioQQQ,
585                                " I could not identify an iso-sequence on this Case A/B command.\n");
586                        fprintf( ioQQQ, " Sorry.\n" );
587                        cdEXIT(EXIT_FAILURE);
588                }
589        }
590        else if( nMatch("DEPA",input.chCARDCAPS ) )
591        {
592                i = 5;
593                /* get expected average departure coefficient, almost certainly 1 */
594                AssertQuantity[nAsserts] =
595                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
596                if( lgEOL )
597                        NoNumb(input.chOrgCard);
598
599                /* this is relative error, max departure from unity of any level or std */
600                AssertError[nAsserts] =
601                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
602                if( lgEOL )
603                        /* default error was set in define above */
604                        AssertError[nAsserts] = ErrorDefault;
605
606                /* H-like key means do one of the hydrogenic ions */
607                if( nMatch("H-LI",input.chCARDCAPS ) )
608                {
609                        /* label is dHII */
610                        strcpy( chAssertLineLabel[nAsserts] , "dHyd" );
611                        /* remember this is departure coefficient for some element */
612                        strcpy( chAssertType[nAsserts] , "D " );
613                        /* now get element number for h ion from element name on card */
614                        if( (wavelength[nAsserts] = (realnum)GetElem(input.chCARDCAPS)) < 0 )
615                        {
616                                fprintf( ioQQQ,
617                                        " I could not identify an element name on this line.\n");
618                                fprintf( ioQQQ, " Sorry.\n" );
619                                cdEXIT(EXIT_FAILURE);
620                        }
621                }
622
623                /* He-like key means do one of the helike ions */
624                else if( nMatch("HE-L",input.chCARDCAPS ) )
625                {
626                        /* label is dHII */
627                        strcpy( chAssertLineLabel[nAsserts] , "d He" );
628                        /* remember this is departure coefficient for some element */
629                        strcpy( chAssertType[nAsserts] , "De" );
630                        /* now get element number for h ion from element name on card */
631                        if( (wavelength[nAsserts] = (realnum)GetElem(input.chCARDCAPS)) < 0 )
632                        {
633                                fprintf( ioQQQ,
634                                        " I could not identify an element name on this line.\n");
635                                fprintf( ioQQQ, " Sorry.\n" );
636                                cdEXIT(EXIT_FAILURE);
637                        }
638                }
639
640                /* this is the large FeII ion */
641                else if( nMatch("FEII",input.chCARDCAPS ) || nMatch("FE II",input.chCARDCAPS ) )
642                {
643                        /* label */
644                        strcpy( chAssertLineLabel[nAsserts] , "d Fe" );
645                        /* remember this is departure coefficient for feii */
646                        strcpy( chAssertType[nAsserts] , "d " );
647                        /* the wavelength is meaningless, but put in 2 since FeII */
648                        wavelength[nAsserts] = 2;
649                }
650
651                /* this is H- h minus */
652                else if( nMatch("HMIN",input.chCARDCAPS ) )
653                {
654                        /* label */
655                        strcpy( chAssertLineLabel[nAsserts] , "d H-" );
656                        /* remember this is departure coefficient for H- */
657                        strcpy( chAssertType[nAsserts] , "d-" );
658                        /* the wavelength is meaningless */
659                        wavelength[nAsserts] = -1;
660                }
661                else
662                {
663                        fprintf( ioQQQ,
664                                " There must be a second key: FEII, H-LIke, HMINus, or HE-Like followed by element.\n");
665                        fprintf( ioQQQ, " Sorry.\n" );
666                        cdEXIT(EXIT_FAILURE);
667                }
668
669                /* last check for key "excited" - which means to skip the ground state */
670                if( nMatch("EXCI",input.chCARDCAPS ) )
671                {
672                        /* this is lowest level - do not do 0 */
673                        AssertQuantity2[nAsserts] = 1.;
674                }
675                else
676                {
677                        /* do the ground state */
678                        AssertQuantity2[nAsserts] = 0.;
679                }
680        }
681
682        /* assert some results from map */
683        else if( nMatch(" MAP",input.chCARDCAPS ) )
684        {
685
686                /* must have heating or cooling, since will check one or the other */
687                /* check heating cooling results from map at some temperature */
688                if( nMatch("HEAT",input.chCARDCAPS ) )
689                {
690                        strcpy( chAssertType[nAsserts] , "mh" );
691                }
692                else if( nMatch("COOL",input.chCARDCAPS ) )
693                {
694                        strcpy( chAssertType[nAsserts] , "mc" );
695                }
696                else
697                {
698                        fprintf( ioQQQ,
699                                " There must be a second key, HEATing or COOLing.\n");
700                        fprintf( ioQQQ, " Sorry.\n" );
701                        cdEXIT(EXIT_FAILURE);
702                }
703
704                /* i was set above for start of scan */
705                /* now get temperature for AssertQuantity2 array*/
706                i = 5;
707                AssertQuantity2[nAsserts] =
708                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
709                if( lgEOL )
710                {
711                        NoNumb(input.chOrgCard);
712                }
713
714                if( AssertQuantity2[nAsserts] <= 10. )
715                {
716                        /* entered as log, but we will compare with linear */
717                        AssertQuantity2[nAsserts] =
718                                pow(10.,AssertQuantity2[nAsserts] );
719                }
720
721                /* print the temperature in the wavelength column */
722                wavelength[nAsserts] = (realnum)AssertQuantity2[nAsserts];
723
724                /* heating or cooling, both log, put into error */
725                AssertQuantity[nAsserts] =
726                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
727                if( lgEOL )
728                {
729                        NoNumb(input.chOrgCard);
730                }
731
732                /* AssertQuantity array will have heating or cooling */
733                AssertQuantity[nAsserts] = pow(10., AssertQuantity[nAsserts]);
734
735                /* optional error, default available (cannot do before loop since we
736                 * do not know how many numbers are on line */
737                AssertError[nAsserts] =
738                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
739                if( lgEOL )
740                {
741                        /* default error was set in define above */
742                        AssertError[nAsserts] = ErrorDefault;
743                }
744
745                /* entered as a log, so print as a log too */
746                lgQuantityLog[nAsserts] = true;
747        }
748
749        /* assert column density of something */
750        else if( nMatch("COLU",input.chCARDCAPS ) )
751        {
752                /* this is column density */
753                strcpy( chAssertType[nAsserts] , "cd" );
754
755                /* this says to look for molecular column density, also could be ion stage */
756                wavelength[nAsserts] = 0;
757
758                /* we want to remember where the match occurred within the string
759                 * since do not want to count the 2 as the first number */
760                /* NB - can't put this expression in the if since many compilers warn against this */
761                if( (i = nMatch(" H2 ",input.chCARDCAPS )) != 0 )
762                {
763                        strcpy( chAssertLineLabel[nAsserts], "H2  " );
764                        /* increment to get past the 2 in the label */
765                        i += 3;
766                        if( nMatch( "LEVE" , input.chCARDCAPS ) )
767                        {
768                                /* this is option for level-specific column density,
769                                 * next two numbers must be v then J */
770                                Param[nAsserts][0] =
771                                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
772                                Param[nAsserts][1] =
773                                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
774                                if( lgEOL )
775                                        NoNumb( input.chCARDCAPS );
776                                /* wavelength will be 10. * vib + rot */
777                                wavelength[nAsserts] = (realnum)(100.*Param[nAsserts][0] + Param[nAsserts][1]);
778                        }
779                        else
780                        {
781                                /* these are flags saying not to do state specific column densities */
782                                Param[nAsserts][0] = -1.;
783                                Param[nAsserts][1] = -1.;
784                        }
785                }
786                else if( (i = nMatch("H3+ ",input.chCARDCAPS )) != 0 )
787                {
788                        strcpy( chAssertLineLabel[nAsserts], "H3+ " );
789                        /* increment to get past the 2 in the label */
790                        i += 3;
791                }
792                else if( (i = nMatch("H2+ ",input.chCARDCAPS )) != 0 )
793                {
794                        strcpy( chAssertLineLabel[nAsserts], "H2+ " );
795                        /* increment to get past the 2 in the label */
796                        i += 3;
797                }
798                else if( (i = nMatch(" H- ",input.chCARDCAPS )) != 0 )
799                {
800                        strcpy( chAssertLineLabel[nAsserts], "H-  " );
801                        /* increment to get past the 2 in the label */
802                        i += 3;
803                }
804                else if( (i = nMatch("H2G ",input.chCARDCAPS )) != 0 )
805                {
806                        strcpy( chAssertLineLabel[nAsserts], "H2g " );
807                        /* increment to get past the 2 in the label */
808                        i += 3;
809                }
810                else if( (i = nMatch("H2* ",input.chCARDCAPS )) != 0 )
811                {
812                        strcpy( chAssertLineLabel[nAsserts], "H2* " );
813                        /* increment to get past the 2 in the label */
814                        i += 3;
815                }
816                else if( (i = nMatch("HEH+",input.chCARDCAPS )) != 0 )
817                {
818                        strcpy( chAssertLineLabel[nAsserts], "HeH+" );
819                        /* increment to get past the 2 in the label */
820                        i += 3;
821                }
822                else if( (i = nMatch(" O2 ",input.chCARDCAPS )) != 0 )
823                {
824                        strcpy( chAssertLineLabel[nAsserts], "O2  " );
825                        /* increment to get past the 2 in the label */
826                        i += 3;
827                }
828                else if( (i = nMatch("H2O ",input.chCARDCAPS )) != 0 )
829                {
830                        strcpy( chAssertLineLabel[nAsserts], "H2O " );
831                        /* increment to get past the 2 in the label */
832                        i += 3;
833                }
834                else if( (i = nMatch(" C2 ",input.chCARDCAPS ) ) !=0 )
835                {
836                        strcpy( chAssertLineLabel[nAsserts], "C2  " );
837                        /* increment to get past the 2 in the label */
838                        i += 3;
839                }
840                else if( (i = nMatch(" C3 ",input.chCARDCAPS ) ) !=0 )
841                {
842                        strcpy( chAssertLineLabel[nAsserts], "C3  " );
843                        /* increment to get past the 3 in the label */
844                        i += 3;
845                }
846                else if( nMatch(" CO ",input.chCARDCAPS ) )
847                {
848                        strcpy( chAssertLineLabel[nAsserts], "CO  " );
849                        i = 5;
850                }
851                else if( nMatch("SIO ",input.chCARDCAPS ) )
852                {
853                        strcpy( chAssertLineLabel[nAsserts], "SiO " );
854                        i = 5;
855                }
856                else if( nMatch(" OH ",input.chCARDCAPS ) )
857                {
858                        strcpy( chAssertLineLabel[nAsserts], "OH  " );
859                        i = 5;
860                }
861                else if( nMatch(" CN ",input.chCARDCAPS ) )
862                {
863                        strcpy( chAssertLineLabel[nAsserts], "CN  " );
864                        i = 5;
865                }
866                else if( nMatch(" CH ",input.chCARDCAPS ) )
867                {
868                        strcpy( chAssertLineLabel[nAsserts], "CH  " );
869                        i = 5;
870                }
871                else if( nMatch(" CH+",input.chCARDCAPS ) )
872                {
873                        strcpy( chAssertLineLabel[nAsserts], "CH+ " );
874                        i = 5;
875                }
876                else
877                {
878                        fprintf( ioQQQ,
879                                " I could not identify H2, H3+, H2+, H2g, H2*, H2H+, CO, O2, SiO, OH, C2 or C3 or on this line.\n");
880                        fprintf( ioQQQ, " Sorry.\n" );
881                        cdEXIT(EXIT_FAILURE);
882                }
883
884                /* i was set above for start of scan */
885                /* now get log of column density */
886                AssertQuantity[nAsserts] =
887                        FFmtRead(input.chCARDCAPS ,&i, INPUT_LINE_LENGTH,&lgEOL);
888                if( lgEOL )
889                {
890                        NoNumb(input.chOrgCard);
891                }
892                /* entered as log, but we will compare with linear */
893                AssertQuantity[nAsserts] =
894                        pow(10.,AssertQuantity[nAsserts] );
895
896                /* optional error, default available (cannot do before loop since we