src/snark/inputfile.cpp

   1 /*
   2  ***********************************************************
   3  $SNARK_Header: S N A R K  1 4 - A PICTURE RECONSTRUCTION PROGRAM $
   4  $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/inputfile.cpp $
   5  $LastChangedRevision: 86 $
   6  $Date: 2014-07-02 16:45:09 -0400 (Wed, 02 Jul 2014) $
   7  $Author: olangthaler $
   8  ***********************************************************
   9
  10  THE FIRST 7 ROWS OF 'OBJECT' STORE THE DENSITY FOR THE 7 ENERGY
  11  LEVELS, 8TH  STORES THE AVERAGE DENSITY;
  12  9TH ROW SPECIFIES TYPE OF OBJECT, 1 FOR ELLIPSE/CIRCLE
  13  2 FOR RECTANGLE/SQUARE
  14  3 FOR TRIANGLE
  15  4 FOR SEGMENT OF CIRCLE
  16  5 FOR SECTOR OF CIRCLE
  17  10TH ROW  CX  X-COORDINATE OF CENTER
  18  11TH ROW  CY  Y-COORDINATE OF CENTER
  19  12TH ROW  U  SEMI X-AXIS LENGTH
  20  13TH ROW  V  SEMI Y-AXIS LENGTH
  21  14TH ROW  ANG  ANGLE OF ROTATION FROM X-AXIS
  22  15TH ROW  SINANG  SIN OF ANGLE OF ROTATION
  23  16TH ROW  COSANG  COSINE OF ANGLE OF ROTATION
  24  17TH ROW RAYL   LENGTH OF RAY INTERSECTED BY OBJECT
  25  THIS VALUE OF SUBSCRIPT REFERRED TO IN CREAPR ONLY
  26  */
  27
  28 #include <cstdlib>
  29 #include <cstdio>
  30 #include <cmath>
  31 #include <time.h>
  32
  33 #include "blkdta.h"
  34 #include "creacm.h"
  35 #include "geom.h"
  36 #include "spctrm.h"
  37 #include "noise.h"
  38 #include "objects.h"
  39 #include "consts.h"
  40 #include "uiod.h"
  41
  42 #include "int2str.h"
  43 #include "anglst.h"
  44 #include "creaer.h"
  45 #include "DIGRand.h"
  46
  47 #include "inputfile.h"
  48
  49 //jk 11/18/2007 introducing variability
  50 REAL variability;
  51 //jk 02/06/09 bug 251
  52 //long newseed;    //the variable Creacm.varseed is used instead
  53 //the array stores values of variability added to each pixel of the phantom
  54 //for each energy level
  55 REAL * variab_array[7];
  56 //the flag is set if the value of variability in the input file is other than zero
  57 BOOLEAN variab_set;
  58 //saves the value of the origical scale which is the mean value for variability calculations
  59 REAL mean;
  60
  61 REAL * phantom_poly_array[7];
  62
  63 void InputFile_class::Open()  //open input file (stdin) and
  64 {
  65         infile = stdin;
  66         if (infile == NULL)
  67         {
  68                 fprintf(output, "\n **** unable to open input file");
  69                 fprintf(output, "\n **** program aborted\n");
  70                 exit(-1);
  71         }
  72 }
  73
  74 void InputFile_class::GetName() // get name
  75 {
  76         int i;
  77         signed CHAR ch;  // changed "CHAR" to "signed CHAR". lajos, Nov 18, 2004
  78
  79         fprintf(output, "\n");
  80
  81         //assign the name of reconstruction to the 'name' member of
  82         //Creacm (defined in global scope in creacm.cpp
  83         for (i = 0; i < 80; i++)
  84         {  // changed "i < 81" to "i < 80". Lajos, Nov 18, 2004
  85                 //if there are fewer than 80 characters, stop earlier
  86                 if (((ch = fgetc(infile)) == '\n') || (ch == EOF))
  87                 {
  88                         break;
  89                 }
  90                 Creacm.name[i] = ch;
  91         }
  92         //skip over remaining characters on the line (only first 80 characters
  93         //are used for the name)
  94         while (ch != '\n' && ch != EOF)
  95                 ch = fgetc(infile);  // bug117 - swr - 6/29/05
  96         // insufficient records in input; EOF encountered
  97         if (ch == EOF)
  98         {
  99                 creaer(25);     //generate an error
 100                 return;         //and quit
 101         }
 102
 103         Creacm.name[i] = 0;     //set the last character to null (it is
 104                                                 //a null terminated string)
 105
 106         fprintf(output, "\n         %s\n", Creacm.name);
 107 }
 108
 109 void InputFile_class::ReadEnergySpectralDistribution()  //READ ENERGY
 110 //SPECTRAL DISTRIBUTION
 111
 112 {
 113         BOOLEAN eol;
 114         INTEGER ndred, tword;
 115         int i;
 116
 117         static const INTEGER spect_codes[3] =
 118         { CHAR2INT('s', 'p', 'e', 'c'), CHAR2INT('m', 'o', 'n', 'o'), CHAR2INT('p',
 119                         'o', 'l', 'y') };
 120
 121         tword = getwrd(TRUE, &eol, spect_codes, 1); /*getwrd returns as its value
 122          the first four letters  of the next string of alpha characters
 123          in the input record  calls function 'getnxt'.
 124          returns eol = TRUE when end of line is reached.*/
 125
 126         if (eol)
 127                 creaer(spect_codes[0]); //generate error if there are
 128         //fewer than 4 characters on the line
 129
 130         Creacm.word2 = getwrd(FALSE, &eol, &(spect_codes[1]), 2); //gets the code
 131         //for either 'mono' or 'poly' from the input file
 132
 133         if (Creacm.word2 == spect_codes[2]) //if 'poly'
 134         {
 135                 // PROCESS POLYCHROMATIC DATA HERE; FOR MONOCHROMATIC GO ELESEWHERE
 136
 137                 Spctrm.nergy = getint(FALSE, &eol);
 138                 if ((Spctrm.nergy <= 0) || (Spctrm.nergy > 7))
 139                         creaer(3);
 140
 141                 getnxt(TRUE);
 142
 143                 for (i = 0; i < Spctrm.nergy; i++)
 144                 {
 145                         Spctrm.energy[i] = getint(FALSE, &eol);
 146                         if (eol)
 147                         {
 148                                 // INSUFFICIENT RECORDS IN INPUT; EOF ENCOUNTERED
 149                                 creaer(25);
 150                                 return;
 151                         };
 152
 153                         Creacm.percnt[i] = getint(FALSE, &eol);
 154                         if (eol)
 155                         {
 156                                 // INSUFFICIENT RECORDS IN INPUT; EOF ENCOUNTERED
 157                                 creaer(25);
 158                                 return;
 159                         }
 160                 };
 161
 162                 fprintf(output,
 163                                 "\n         energy spectrum is polychromatic with %3i     energy levels",
 164                                 Spctrm.nergy);
 165                 fprintf(output, "\n         with the following photon distribution");
 166                 fprintf(output, "\n         energy level   percent");
 167                 for (i = 0; i < Spctrm.nergy; i++)
 168                 {
 169                         fprintf(output, "\n%17i %11i", Spctrm.energy[i], Creacm.percnt[i]);
 170                 }
 171
 172                 // COMPUTE THE WEIGHTS FOR EACH ENERGY LEVEL FOR AVERAGING THE
 173                 // ABSORPTION OVER THE ENERGY LEVELS
 174                 ndred = 0;
 175
 176                 for (i = 0; i < Spctrm.nergy; i++)
 177                 {
 178                         Spctrm.engwt[i] = Creacm.percnt[i] * (REAL) 0.01;
 179                         ndred += Creacm.percnt[i];
 180                         if (Creacm.percnt[i] < 0)
 181                                 ndred = -100000;
 182                 }
 183
 184                 if (ndred != 100)
 185                         creaer(4);
 186         }
 187         else    //if 'mono'
 188         {
 189                 if (Creacm.word2 == spect_codes[1])
 190                 {
 191                         // MONOCHROMATIC DATA
 192
 193                         Spctrm.nergy = 1;
 194                         Creacm.percnt[0] = 100;
 195                         Spctrm.engwt[0] = 1.0;
 196                         Spctrm.energy[0] = getint(FALSE, &eol);
 197                         ;
 198                         fprintf(output,
 199                                         "\n         energy spectrum is monochromatic at energy level %5i\n",
 200                                         Spctrm.energy[0]);
 201                 }
 202                 else
 203                         creaer(21);
 204         }
 205 }
 206
 207 void InputFile_class::GetInformationABoutCreationOfPhantom(INTEGER invoke)
 208 // GET INFORMATION ABOUT CREATION OF PHANTOM
 209 {
 210         BOOLEAN eol;
 211         static const INTEGER crea = CHAR2INT('c', 'r', 'e', 'a');
 212
 213         INTEGER tword;
 214         static const INTEGER phan_codes[2] =
 215         { CHAR2INT('p', 'h', 'a', 'n'), CHAR2INT('a', 'v', 'e', 'r') };
 216
 217         tword = getwrd(TRUE, &eol, phan_codes, 1);
 218         if (eol)
 219                 creaer(phan_codes[0]);
 220
 221         fprintf(output, "\n");
 222
 223         Creacm.word2 = getwrd(FALSE, &eol, &(phan_codes[1]), 1);
 224         Creacm.picflg = (Creacm.word2 == phan_codes[1]);
 225
 226         // IF PHANTOM NOT CREATED DONT READ NAVE1
 227
 228         if (Creacm.picflg)
 229         {
 230                 if (invoke == crea)
 231                 {
 232                         fprintf(output, "\n         this run will generate a phantom");
 233                 }
 234                 Creacm.nave1 = getint(FALSE, &eol);
 235                 fprintf(output,
 236                                 "\n         density in each pixel is obtained as the average of %i x %i points\n",
 237                                 Creacm.nave1, Creacm.nave1);
 238
 239                 if ((Creacm.nave1 <= 0) || (Creacm.nave1 == Creacm.nave1 / 2 * 2))
 240                         creaer(8);
 241
 242                 // READ IN GRID SIZE AND PIXEL SIZE .
 243
 244                 GeoPar.nelem = getint(TRUE, &eol);
 245                 GeoPar.pixsiz = getnum(FALSE, &eol);
 246                 fprintf(output, "\n         picture size %i x %i,  pixel size %10.4f\n",
 247                                 GeoPar.nelem, GeoPar.nelem, GeoPar.pixsiz);
 248
 249                 if ((GeoPar.nelem <= 0) || (GeoPar.nelem == GeoPar.nelem / 2 * 2))
 250                         creaer(9);
 251                 if (GeoPar.pixsiz < Consts.zero)
 252                         creaer(10);
 253                 GeoPar.midpix = (GeoPar.nelem + 1) / 2;
 254                 GeoPar.area = GeoPar.nelem * GeoPar.nelem;
 255         }
 256
 257 }
 258
 259 // READ BACKGROUND ABSORPTION
 260 void InputFile_class::ReadBackgroundAbsorption()
 261 {
 262         BOOLEAN eol;
 263         int i;
 264         static const INTEGER hback = CHAR2INT('b', 'a', 'c', 'k');
 265
 266         if (word != hback)
 267                 creaer(hback);
 268         for (i = 0; i < Spctrm.nergy; i++)
 269         {
 270                 Spctrm.backgr[i] = getnum(FALSE, &eol);
 271         }
 272
 273         fprintf(output, "\n                               at levels");
 274         fprintf(output, "\n                               ");
 275         for (i = 0; i < Spctrm.nergy; i++)
 276         {
 277                 fprintf(output, "%9i ", Spctrm.energy[i]);
 278         }
 279
 280         fprintf(output, "\n          background absorption");
 281         for (i = 0; i < Spctrm.nergy; i++)
 282         {
 283                 fprintf(output, "%9.4f ", Spctrm.backgr[i]);
 284         }
 285 }
 286
 287 // READ RAYSUM INFORMATION
 288 void InputFile_class::ReadRaysumInformation()
 289 {
 290         BOOLEAN eol;
 291         //  static const INTEGER haver = CHAR2INT('a','v','e','r');  // commented out since it is unused. Lajos, Jan 13, 2005
 292
 293         INTEGER tword;
 294         static const INTEGER rays_codes[2] =
 295         { CHAR2INT('r', 'a', 'y', 's'), CHAR2INT('a', 'v', 'e', 'r') };
 296
 297         tword = getwrd(TRUE, &eol, rays_codes, 1);
 298         if (eol)
 299                 creaer(rays_codes[0]);
 300
 301         fprintf(output, "\n");
 302
 303         Creacm.word2 = getwrd(FALSE, &eol, &(rays_codes[1]), 1);
 304         Creacm.prjflg = (Creacm.word2 == rays_codes[1]);
 305 }
 306
 307 // INPUT ENTIRE GEOMETRY DATA
 308 void InputFile_class::ReadGeometryData(INTEGER invoke)
 309 {
 310         BOOLEAN eol;
 311         INTEGER tword;
 312
 313         static const INTEGER crea = CHAR2INT('c', 'r', 'e', 'a');
 314
 315         static const INTEGER geom_codes[10] =  /// bug #135. Lajos, 08/01/2005
 316                         {
 317                                         CHAR2INT('g', 'e', 'o', 'm'), CHAR2INT('p', 'a', 'r', 'a'),
 318                                         CHAR2INT('d', 'i', 'v', 'e'), CHAR2INT('l', 'i', 'n', 'o'), /// bug #135. Lajos, 08/01/2005
 319                                         CHAR2INT('u', 'n', 'i', 'f'), CHAR2INT('v', 'a', 'r', 'i'),
 320                                         CHAR2INT('l', 'i', 'n', 'e'), CHAR2INT('s', 't', 'r', 'i'),
 321                                         CHAR2INT('a', 'r', 'c', ' '), CHAR2INT('t', 'a', 'n', 'g') };
 322
 323         fprintf(output, "\n");
 324
 325         tword = getwrd(TRUE, &eol, geom_codes, 1);
 326         if (eol)
 327                 creaer(geom_codes[0]);
 328
 329         fprintf(output, "\n");
 330
 331         // GET RAY TYPE, PARALLEL, DIVERGENT, OR LINOGRAM, UNIF, VARI, TANG OR ARC  /// bug #135. Lajos, 08/01/2005
 332         tword = getwrd(TRUE, &eol, &(geom_codes[1]), 3); /// bug #135. Lajos, 08/01/2005
 333         if (eol)
 334         {
 335                 fprintf(output,
 336                                 "\n **** keyword parallel, divergent, or linogram is missing"); /// bug #135. Lajos, 08/01/2005
 337                 fprintf(output, "\n **** program aborted\n");
 338                 exit(-1);
 339         };
 340         Creacm.word1 = tword;
 341
 342         GeoPar.lino = FALSE;  /// bug #135. Lajos, 08/02/2005
 343
 344         if (tword == geom_codes[1])
 345         {
 346                 tword = getwrd(FALSE, &eol, &(geom_codes[4]), 2); /// bug #135. Lajos, 08/01/2005
 347                 if (eol)
 348                 {
 349                         fprintf(output, "\n **** keyword uniform or variable is missing");
 350                         fprintf(output, "\n **** program aborted\n");
 351                         exit(-1);
 352                 }
 353         }
 354         else if (tword == geom_codes[2])  /// bug #135. Lajos, 08/01/2005
 355         {
 356                 tword = getwrd(FALSE, &eol, &(geom_codes[8]), 2); /// bug #135. Lajos, 08/01/2005
 357                 if (eol)
 358                 {
 359                         fprintf(output, "\n **** keyword arc or tangent is missing");
 360                         fprintf(output, "\n **** program aborted\n");
 361                         exit(-1);
 362                 }
 363         }
 364         else  /// bug #135. Lajos, 08/01/2005
 365         {
 366                 // assuming PARALLEL VARIABLE LINE when LINOGRAM has been specified
 367                 GeoPar.lino = TRUE;
 368                 Creacm.word1 = geom_codes[1];
 369                 tword = geom_codes[5];
 370         }
 371
 372         Creacm.word2 = tword;
 373
 374         GeoPar.par = (Creacm.word1 == geom_codes[1]);
 375         GeoPar.div = (Creacm.word1 == geom_codes[2]); // moved out of the "if(!GeoPar.par) {...}" block. Lajos, 08/02/2005
 376
 377         if (!GeoPar.par)
 378         {
 379
 380                 // DIVERGENT GEOMETRY READ SPECS
 381
 382                 GeoPar.line = TRUE;
 383                 GeoPar.arc = (Creacm.word2 == geom_codes[8]); /// bug #135. Lajos, 08/01/2005
 384                 GeoPar.tang = (Creacm.word2 == geom_codes[9]); /// bug #135. Lajos, 08/01/2005
 385                 if (!(GeoPar.div && (GeoPar.arc || GeoPar.tang)))
 386                         creaer(11);
 387
 388                 fprintf(output, "\n         rays are divergent from point sources");
 389
 390                 GeoPar.radius = getnum(FALSE, &eol);
 391                 GeoPar.stod = getnum(FALSE, &eol);
 392                 fprintf(output, "\n         source to origin distance  %10.4f",
 393                                 GeoPar.radius);
 394
 395                 if (GeoPar.arc)
 396                 {
 397                         fprintf(output,
 398                                         "\n         the detectors lie on an arc with source to detector distance = %10.4f",
 399                                         GeoPar.stod);
 400                 }
 401
 402                 if (GeoPar.tang)
 403                 {
 404                         fprintf(output,
 405                                         "\n         the detectors lie on a st. line with source to middle detector distance = %10.4f",
 406                                         GeoPar.stod);
 407                 }
 408
 409                 if (invoke != crea)
 410                 {
 411                         if (GeoPar.radius <= GeoPar.nelem * GeoPar.pixsiz / 1.414)
 412                                 creaer(11);
 413                 }
 414         }
 415         else
 416         {
 417                 // PARALLEL GEOMETRY  READ SPECS; IS SPACING UNIFORM OR VARIABLE
 418
 419                 GeoPar.uni = (Creacm.word2 == geom_codes[4]); /// bug #135. Lajos, 08/01/2005
 420                 GeoPar.vri = (Creacm.word2 == geom_codes[5]); /// bug #135. Lajos, 08/01/2005
 421
 422                 if (GeoPar.lino)
 423                         tword = geom_codes[6];  /// bug #135. Lajos, 08/01/2005
 424                 else
 425                 {
 426                         tword = getwrd(FALSE, &eol, &(geom_codes[6]), 2); /// bug #135. Lajos, 08/01/2005
 427                         if (eol)
 428                         {
 429                                 fprintf(output, "\n **** keyword line or strip is missing");
 430                                 fprintf(output, "\n **** program aborted\n");
 431                                 exit(-1);
 432                         }
 433                 }
 434                 Creacm.word3 = tword;
 435
 436                 GeoPar.line = (Creacm.word3 == geom_codes[6]); /// bug #135. Lajos, 08/01/2005
 437                 GeoPar.strip = (Creacm.word3 == geom_codes[7]); /// bug #135. Lajos, 08/01/2005
 438
 439                 if (!(GeoPar.par && (GeoPar.strip || GeoPar.line)))
 440                         creaer(11);
 441                 if (!(GeoPar.uni || GeoPar.vri))
 442                         creaer(12);
 443
 444                 fprintf(output, "\n         rays are parallel with");
 445                 if (GeoPar.uni)
 446                 {
 447                         fprintf(output, " uniform spacing between rays");
 448                 }
 449                 if (GeoPar.vri)
 450                 {
 451                         fprintf(output, " variable spacing between rays");
 452                 }
 453                 if (GeoPar.strip)
 454                 {
 455                         fprintf(output, "\n         data collected along strips\n");
 456                 }
 457                 if (GeoPar.line)
 458                 {
 459                         fprintf(output, "\n         data collected along lines\n");
 460                 }
 461         }
 462         fprintf(output, "\n");
 463 }
 464
 465 void InputFile_class::ReadNumberOfRaysAndRaySpacing(INTEGER invoke)
 466 {
 467         BOOLEAN eol;
 468         INTEGER melen, tword;
 469         REAL zisxip;
 470
 471         static const INTEGER proj = CHAR2INT('p', 'r', 'o', 'j');
 472         static const INTEGER crea = CHAR2INT('c', 'r', 'e', 'a');
 473
 474         static const INTEGER rays_codes[3] =
 475         { CHAR2INT('r', 'a', 'y', 's'), CHAR2INT('u', 's', 'e', 'r'), CHAR2INT('p',
 476                         'r', 'o', 'g') };
 477
 478         if (GeoPar.lino)  /// bug #135. Lajos, 08/01/2005
 479         {
 480                 // assuming RAYS USER if LINOGRAM has been specified
 481                 if (GeoPar.nelem == 0)
 482                 {
 483                         fprintf(output,
 484                                         "\n **** nelem must be specified using phantom average if linogram geometry has been chosen");
 485                         fprintf(output, "\n **** program aborted\n");
 486                         exit(-1);
 487                 }
 488                 tword = rays_codes[1];
 489                 GeoPar.usrays = 2 * GeoPar.nelem + 1;
 490                 GeoPar.pinc = GeoPar.pixsiz;
 491         }
 492         else
 493         {
 494                 tword = getwrd(TRUE, &eol, rays_codes, 1);
 495                 if (eol)
 496                         creaer(rays_codes[0]);
 497
 498                 GeoPar.usrays = 0;
 499
 500                 tword = getwrd(FALSE, &eol, &(rays_codes[1]), 2);
 501                 if (eol)
 502                         creaer(15);
 503
 504                 if (tword == rays_codes[1])
 505                 {
 506
 507                         GeoPar.usrays = getint(FALSE, &eol);
 508
 509                         if (!((GeoPar.usrays == 0) && (invoke == crea)))
 510                         {
 511                                 if ((GeoPar.usrays % 2) <= 0)
 512                                         creaer(13);
 513                         }
 514                 }
 515                 else
 516                 {
 517                         melen = getint(FALSE, &eol);
 518                         zisxip = getnum(FALSE, &eol);
 519                 }
 520
 521                 GeoPar.pinc = getnum(FALSE, &eol);
 522                 if (GeoPar.pinc < Consts.zero)
 523                         creaer(10);
 524         }
 525
 526         if ((invoke == proj) && (GeoPar.nelem == 0))
 527                 creaer(23);
 528         if (GeoPar.arc)
 529         {
 530                 if ((GeoPar.pinc * GeoPar.usrays - 1.0) / GeoPar.stod > Consts.pi)
 531                         creaer(11);
 532         }
 533
 534         if (tword == rays_codes[2])
 535                 GeoPar.usrays = GeoPar.numray(melen, zisxip);
 536
 537         if (GeoPar.nelem > 0)
 538                 GeoPar.snrays = GeoPar.numray(GeoPar.nelem, GeoPar.pixsiz);
 539         if (GeoPar.nelem == 0)
 540                 GeoPar.snrays = GeoPar.usrays;
 541         if (GeoPar.usrays == 0)
 542                 GeoPar.usrays = GeoPar.snrays;
 543
 544         fprintf(output, "\n         number of rays per projection %5i",
 545                         GeoPar.usrays);
 546
 547         if (invoke == proj)
 548         {
 549                 fprintf(output, "\n         snark computed number of rays %5i",
 550                                 GeoPar.snrays);
 551         }
 552         fprintf(output, "\n         at detector spacing %10.4f", GeoPar.pinc);
 553
 554         // SET THE APPROPRIATE VALUES FOR THE FIRST AND LAST RAY NUMBERS
 555         // FOR THE USER RAY SET AND THE SNARK RAY SET
 556         // THE VALUE OF NRAYS IN /USER/ COMMON NOW IS NO OF USER SPECIFIED
 557         // RAYS IF GIVEN  OR  SNRAYS COMPUTED BY SNARK IF DEFAULT
 558
 559         GeoPar.nrays = MAX0(GeoPar.snrays, GeoPar.usrays);
 560         GeoPar.midray = GeoPar.nrays / 2;
 561
 562         GeoPar.fsnray = (GeoPar.nrays - GeoPar.snrays) / 2;
 563         GeoPar.lsnray = GeoPar.fsnray + GeoPar.snrays - 1;
 564
 565         GeoPar.fusray = (GeoPar.nrays - GeoPar.usrays) / 2;
 566         GeoPar.lusray = GeoPar.fusray + GeoPar.usrays - 1;
 567 }
 568
 569 // MEASUREMENT STATISTICS  --  NOISE AND SCATTER
 570
 571 void InputFile_class::GetMesurementStatistics()
 572 {
 573         BOOLEAN eol;
 574         BOOLEAN perfct, noisfl;
 575
 576         int n;
 577
 578         static const INTEGER hquan = CHAR2INT('q', 'u', 'a', 'n');
 579         static const INTEGER hscat = CHAR2INT('s', 'c', 'a', 't');
 580         static const INTEGER hmult = CHAR2INT('m', 'u', 'l', 't');
 581         static const INTEGER haddi = CHAR2INT('a', 'd', 'd', 'i');
 582
 583         static const INTEGER meas_codes[10] =
 584         {
 585                         CHAR2INT('m', 'e', 'a', 's'), CHAR2INT('p', 'e', 'r', 'f'),
 586                         CHAR2INT('n', 'o', 'i', 's'), CHAR2INT('s', 'e', 'e', 'd'),
 587                         CHAR2INT('b', 'a', 'c', 'k'), CHAR2INT('q', 'u', 'a', 'n'),
 588                         CHAR2INT('s', 'c', 'a', 't'), CHAR2INT('m', 'u', 'l', 't'),
 589                         CHAR2INT('a', 'd', 'd', 'i'), CHAR2INT('c', 'a', 'l', 'i')
 590         };
 591
 592         NoisePar.setQuantumNoise(0, 1.0, 1.0);  //bug221 - swr - 03/03/07
 593
 594         word = getwrd(TRUE, &eol, meas_codes, 1);
 595         if (eol)
 596                 creaer(meas_codes[0]);
 597
 598         word = getwrd(FALSE, &eol, &(meas_codes[1]), 2);
 599         if (eol)
 600         {
 601                 fprintf(output, "\n **** keyword perfect or noisy is missing");
 602                 fprintf(output, "\n **** program aborted\n");
 603                 exit(-1);
 604         };
 605
 606         perfct = (word == meas_codes[1]);
 607         noisfl = (word == meas_codes[2]);
 608         if (!noisfl)
 609         {
 610
 611                 // IF PERFECT MEASUREMENTS THEN SKIP NOISE INFORMATION
 612                 if (!perfct)
 613                         creaer(meas_codes[1]);
 614
 615                 fprintf(output, "\n         projection data are noiseless\n");
 616
 617                 word = getwrd(TRUE, &eol, &(meas_codes[4]), 1);
 618         }
 619         else
 620         {
 621                 // NOISE INFORMATION
 622
 623                 fprintf(output,
 624                                 "\n         noise characteristics of projection data follow");
 625                 fprintf(output, "\n               nature          characteristics");
 626
 627                 for (;;)
 628                 {
 629                         // CHECK TYPE OF NOISE AND READ RELEVANT DATA
 630                         word = getwrd(TRUE, &eol, &(meas_codes[3]), 6);
 631
 632                         if ((word == meas_codes[3]) || (word == meas_codes[4]))
 633                                 break;
 634
 635                         switch (word)
 636                         {
 637
 638                         case hquan:
 639                                 // QUANTUM NOISE
 640                                 NoisePar.quanmn = getnum(FALSE, &eol);
 641                                 NoisePar.quancm = getnum(FALSE, &eol);
 642
 643                                 // GET CALIBRATION TYPE AND ECHO ON OUTPUT
 644                                 NoisePar.quanin = 1;
 645                                 if (getwrd(FALSE, &eol, &(meas_codes[9]), 1) == meas_codes[9])
 646                                         NoisePar.quanin = getint(FALSE, &eol);
 647                                 if (NoisePar.quanin == 4)
 648                                 {
 649                                         fprintf(output, "\n              Emission tomography");
 650                                 }
 651                                 else
 652                                 {
 653
 654                                         if ((NoisePar.quanin < 1) || (NoisePar.quanin > 4))
 655                                                 creaer(16);
 656                                         fprintf(output,
 657                                                         "\n               quantum        mean number of photons in actual measurement %12.4e",
 658                                                         NoisePar.quanmn);
 659                                         fprintf(output,
 660                                                         "\n                              mean number of photons in calibration measurement is %10.4f    times above",
 661                                                         NoisePar.quancm);
 662
 663                                         if (NoisePar.quanin == 1)
 664                                         {
 665                                                 fprintf(output,
 666                                                                 "\n                              calibration is constant for all rays in any single projection");
 667                                         }
 668
 669                                         if (NoisePar.quanin == 2)
 670                                         {
 671                                                 fprintf(output,
 672                                                                 "\n                              calibration is constant over all projections for any single ray");
 673                                         }
 674
 675                                         if (NoisePar.quanin == 3)
 676                                         {
 677                                                 fprintf(output,
 678                                                                 "\n                              calibration is variable for every ray in every projection");
 679                                         }
 680
 681                                         if ((NoisePar.quanmn < Consts.zero)
 682                                                         || (NoisePar.quancm < Consts.zero))
 683                                                 creaer(26);
 684
 685                                         // DISTRIBUTE PHOTONS ACCORDING TO APERTURE WEIGHTS COMPUTED EARLIER
 686
 687                                         for (n = 0; n < GeoPar.nave2; n++)
 688                                         {
 689                                                 Creacm.aperwt[n] *= NoisePar.quanmn;
 690                                         }
 691
 692                                         fprintf(output,
 693                                                         "\n                              mean numbers of photons entering each substrip are");
 694                                         fprintf(output, "\n                              ");
 695                                         for (n = 0; n < GeoPar.nave2; n++)
 696                                         {
 697                                                 fprintf(output, "%12.4e ", Creacm.aperwt[n]);
 698                                         }
 699                                 }
 700                                 NoisePar.setQuantumNoise();  //bug221 - swr - 03/03/07
 701                                 //}
 702                                 break;
 703
 704                         case hscat:
 705                                 ///if(word == hscat) {
 706                                 // SCATTER NOISE
 707                                 //bug221 - swr - 03/03/07
 708                         {
 709                                 REAL peak = getnum(FALSE, &eol);
 710                                 REAL width = getnum(FALSE, &eol);
 711                                 NoisePar.setScatter(peak, width);
 712                                 fprintf(output,
 713                                                 "\n               scatter        the fraction scattered in line is %10.4f",
 714                                                 peak);
 715                                 fprintf(output,
 716                                                 "\n                              scatter extends to a distance of %10.4f",
 717                                                 width);
 718                                 fprintf(output,
 719                                                 "\n                              dying linearly from the middle");
 720
 721                                 if ((peak < Consts.zero) || (width < Consts.zero))
 722                                         creaer(17);
 723                                 break;
 724                         }
 725                                 ///}
 726
 727                         case hmult:
 728                                 // MULTIPLICATIVE NOISE
 729                                 //bug221 - swr - 03/03/07
 730                         {
 731                                 REAL mean = getnum(FALSE, &eol);
 732                                 REAL std = getnum(FALSE, &eol);
 733                                 NoisePar.setMultiplicativeNoise(mean, std);
 734                                 fprintf(output,
 735                                                 "\n               multiplicative mean %10.4f     std. dev. %10.4f",
 736                                                 mean, std);
 737                                 if (fabs(NoisePar.ultnmn) < Consts.zero)
 738                                         creaer(27);
 739                                 break;
 740                         }
 741
 742                         case haddi:
 743                                 // ADDITIVE NOISE
 744                                 //bug221 - swr - 03/03/07
 745                         {
 746                                 REAL mean = getnum(FALSE, &eol);
 747                                 REAL std = getnum(FALSE, &eol);
 748                                 NoisePar.setAdditiveNoise(mean, std);
 749                                 fprintf(output,
 750                                                 "               additive       %10.4f     std. dev. %10.4f",
 751                                                 mean, std);
 752                                 break;
 753                         }
 754
 755                         default:
 756                                 creaer(18);
 757                         }
 758                 }
 759                 noisfl = ((NoisePar.quanin > 0) || NoisePar.ultnfl || NoisePar.addnfl);
 760
 761                 if (noisfl)
 762                 {
 763                         if (word != meas_codes[3])
 764                                 creaer(meas_codes[3]);
 765
 766                         // INPUT SEED FOR RANDOM NUMBER GENERATOR
 767
 768                         Creacm.iseed = getint(FALSE, &eol);
 769                         fprintf(output,
 770                                         "\n          seed for random number generator is     %10i\n",
 771                                         Creacm.iseed);
 772
 773                         word = getwrd(TRUE, &eol, &(meas_codes[4]), 1);
 774                 }
 775         }
 776 }
 777
 778 // READ NUMBER OF PROJECTIONS
 779 // READ LIST OF PROJECTION ANGLES INTO WORK AREA
 780
 781 void InputFile_class::ReadProjectionAngles(BOOLEAN f11)
 782 {
 783         BOOLEAN eol;
 784         INTEGER JJ, IANG, NUMRAYS, II, LIMIT, NN;
 785         REAL RAD_DEG;
 786         int i, j, np;
 787         REAL delta;
 788         INTEGER tword;
 789
 790         BOOLEAN equisp;
 791         BOOLEAN uequisp;
 792         BOOLEAN linosp;
 793
 794         static const INTEGER angl_codes[2] =  /// bug #135. Lajos, 08/01/2005
 795                         { CHAR2INT('a', 'n', 'g', 'l'), CHAR2INT('e', 'q', 'u', 'a') /// bug #135. Lajos, 08/01/2005
 796                         //    CHAR2INT('l','i','n','o')  /// bug #135. Lajos, 08/01/2005
 797                         };
 798
 799         fprintf(output, "\n");
 800
 801         if (GeoPar.lino)  /// bug #135. Lajos, 08/01/2005
 802         {
 803                 GeoPar.prjnum = 2 * (2 * GeoPar.nelem + 1);
 804         }
 805         else
 806         {
 807                 tword = getwrd(TRUE, &eol, angl_codes, 1);
 808                 if (eol)
 809                         creaer(angl_codes[0]);
 810
 811                 GeoPar.prjnum = getint(FALSE, &eol);
 812         }
 813         fprintf(output, "\n         total number of projections %5i",
 814                         GeoPar.prjnum);
 815
 816         if (GeoPar.prjnum <= 0)
 817                 creaer(14);
 818
 819         // FIRST ALLOCATE STORAGE AREA
 820         if (Creacm.pang != NULL)
 821                 delete[] Creacm.pang;
 822         Creacm.pang = new REAL[GeoPar.prjnum];
 823         Anglst.Init(GeoPar.prjnum);
 824
 825         // WANT TO GENERATE EQUALLY-SPACED ANGLES ONLY IF SPECIFIED
 826
 827         equisp = FALSE;
 828         linosp = FALSE;
 829         uequisp = FALSE;
 830
 831         if (GeoPar.lino)
 832                 linosp = TRUE;  /// bug #135. Lajos, 08/01/2005
 833         else
 834         {
 835                 tword = getwrd(FALSE, &eol, &(angl_codes[1]), 1); /// bug #135. Lajos, 08/01/2005
 836                 if (tword == angl_codes[1])
 837                 {
 838                         equisp = TRUE;
 839                 }
 840                 else
 841                 {
 842                         uequisp = TRUE;
 843                 }
 844         }
 845
 846         if (equisp)
 847         {
 848                 getnxt(TRUE);
 849
 850                 *Creacm.pang = getnum(FALSE, &eol);
 851                 if (eol)
 852                         creaer(25);  // added error checking. Lajos, 03/29/2005
 853                 *(Creacm.pang + 1) = getnum(FALSE, &eol);
 854                 if (eol)
 855                         creaer(25);  // added error checking. Lajos, 03/29/2005
 856                 //C     intang = intang + 1
 857
 858                 // CHANGED CODE FOR READING EQUALLY SPACED ANGLES TO MAKE USERS
 859                 // LIFE BIT MORE COMFORTABLE FOR VERSION 1.02 SNARK87 10th AUG
 860                 // FOR flag DESCRIPTION SEE CODE AT LABEL 40
 861                 // -----S.SADANANDA AUG 10 1987
 862
 863                 if (GeoPar.prjnum <= 1)
 864                         creaer(20);
 865
 866                 delta = (Creacm.pang[1] - Creacm.pang[0]) / (GeoPar.prjnum - 1);
 867 #ifdef FFCOMPARE
 868                 fprintf(output,"\n %36.30f", (double) delta);
 869 #endif
 870                 for (np = 1; np < GeoPar.prjnum; np++)
 871                 {
 872                         Creacm.pang[np] = np * delta + Creacm.pang[0];
 873 #ifdef FFCOMPARE
 874                         fprintf(output,"\n %5i %36.30f", np+1, Creacm.pang[np]);
 875 #endif
 876                 }
 877         }
 878
 879         // Modified code to enable generation of angles for linogram
 880         // reconstruction
 881         // ----- J. Browne, May 1993
 882
 883         else if (linosp)
 884         {
 885                 if (!(GeoPar.par && GeoPar.vri))
 886                 {
 887                         fprintf(output,
 888                                         "\n **** for linogram reconstruction use PARALLEL geometry with VARIABLE detector spacing");
 889                         fprintf(output, "\n **** program aborted\n");
 890                         exit(-1);
 891                 }
 892                 NN = (GeoPar.prjnum / 2 - 3) / 4;
 893                 LIMIT = 2 * NN + 1;
 894                 NUMRAYS = 4 * NN + 3;
 895                 RAD_DEG = (REAL) 180.0 / Consts.pi;
 896                 JJ = 0;
 897                 for (II = 1; II <= 2; II++)
 898                 {
 899                         for (IANG = -LIMIT; IANG <= LIMIT; IANG++)
 900                         {
 901
 902                                 if (II == 1)
 903                                 {
 904                                         Creacm.pang[JJ] = RAD_DEG
 905                                                         * (REAL) atan((REAL) 2.0 * IANG / NUMRAYS) + 90;
 906                                 }
 907
 908                                 if (II == 2)
 909                                 {
 910                                         Creacm.pang[JJ] = RAD_DEG
 911                                                         * (REAL) atan((REAL) 2.0 * IANG / NUMRAYS) + 180;
 912                                 }
 913                                 JJ++;
 914
 915                         }
 916                 }
 917         }
 918         else if (uequisp)
 919         {
 920                 if (f11)
 921                 { // from file11
 922
 923                         eol = FALSE;
 924                         getnxt(TRUE);
 925                         fprintf(output, "\n");
 926                         //bug 276, jklukowska, allow free format in file11
 927                         for (i = 0; i < GeoPar.prjnum; i++)
 928                         {
 929                                 Creacm.pang[i] = getnum(FALSE, &eol);    // try to get number
 930                                 while (eol)
 931                                 {                    // if end of line
 932                                         getnxt(TRUE, FALSE);         // get another line
 933                                         Creacm.pang[i] = getnum(FALSE, &eol);  // try again
 934                                 }
 935                         }
 936
 937                 }
 938                 else
 939                 {  // from input file
 940                         eol = FALSE;
 941                         getnxt(TRUE);
 942
 943                         for (i = 0; i < GeoPar.prjnum; i++)
 944                         {
 945                                 Creacm.pang[i] = getnum(FALSE, &eol);    // try to get number
 946                                 while (eol)
 947                                 {                    // if end of line
 948                                         getnxt(TRUE);         // get another line
 949                                         Creacm.pang[i] = getnum(FALSE, &eol);  // try again
 950                                 }
 951                         }
 952                 }
 953         }
 954
 955         fprintf(output, "\n         projection angles");
 956
 957         for (j = 0; j < GeoPar.prjnum; j++)
 958         {
 959                 if (((j % 10) == 0) && (j > 0))
 960                         fprintf(output, "\n                          ");
 961                 fprintf(output, " %9.4f", Creacm.pang[j]);
 962         }
 963
 964         // CONVERT TO RADIANS AND COMPUTE SIN AND COS TABLES
 965
 966         Anglst.InitDiv(Creacm.pang, GeoPar.prjnum);
 967
 968         fprintf(output, "\n");
 969 }
 970
 971 void InputFile_class::ReadIfoForComputionRaysums()
 972 {
 973         BOOLEAN eol;
 974         int i, n, ns;
 975
 976         GeoPar.nave2 = getint(FALSE, &eol);
 977         fprintf(output,
 978                         "\n         projection data are calculated by dividing each ray interval into %i substrips",
 979                         GeoPar.nave2);
 980
 981         if ((GeoPar.nave2 <= 0) || (GeoPar.nave2 > 13)
 982                         || (GeoPar.nave2 == GeoPar.nave2 / 2 * 2))
 983                 creaer(8);
 984         getnxt(TRUE);
 985
 986         /////////////////////////////////////////////////////////////
 987         for (i = 0; i < GeoPar.nave2; i++)
 988         { /// Is this what is needed MK ???
 989                 GeoPar.naper[i] = getint(FALSE, &eol);
 990                 if (eol == TRUE)
 991                 {
 992                         // INSUFFICIENT RECORDS IN INPUT; EOF ENCOUNTERED
 993                         creaer(25);
 994                 }
 995         }
 996
 997         fprintf(output, "\n         with aperture (substrip) weights");
 998         for (n = 0; n < GeoPar.nave2; n++)
 999         {
1000                 fprintf(output, " %5i", GeoPar.naper[n]);
1001         }
1002
1003         // GET SUM OF THE APERTURE WEIGHTS FOR LATER USE
1004         ns = 0;
1005         for (n = 0; n < GeoPar.nave2; n++)
1006         {
1007                 if (GeoPar.naper[n] < 0)
1008                         creaer(19);
1009                 ns += GeoPar.naper[n];
1010         }
1011         if (ns <= 0)
1012                 creaer(19);
1013
1014         // COMPUTE NORMALIZED APERTURE WEIGHTS FOR USE IN CREATING
1015         // PROJECTION DATA
1016         for (n = 0; n < GeoPar.nave2; n++)
1017         {
1018                 Creacm.aperwt[n] = ((REAL) GeoPar.naper[n]) / ((REAL) (ns));
1019         }
1020 }
1021
1022 // READ OBJECTS MAKING UP THE PHANTOM
1023
1024 void InputFile_class::ReadObjects(INTEGER invoke)
1025 {
1026         int otype;
1027         int i, j, n;
1028
1029         REAL densit;
1030
1031         INTEGER itype, tword;
1032         REAL cxt;
1033         REAL cyt;
1034         REAL ut;
1035         REAL vt;
1036         REAL angt;
1037
1038         REAL scale;
1039         REAL phi;
1040         BOOLEAN eol;
1041         REAL dent[7];
1042
1043         SNARK_Object_Node *f_object, *t_object;
1044
1045         static const INTEGER crea = CHAR2INT('c', 'r', 'e', 'a');
1046
1047         static const INTEGER object_codes[8] =
1048         {
1049                         CHAR2INT('o', 'b', 'j', 'e'), CHAR2INT('e', 'l', 'i', 'p'),
1050                         CHAR2INT('r', 'e', 'c', 't'), CHAR2INT('t', 'r', 'i', 'a'),
1051                         CHAR2INT('s', 'e', 'g', 'm'), CHAR2INT('s', 'e', 'c', 't'),
1052                         CHAR2INT('l', 'a', 's', 't'), CHAR2INT('d', 'e', 'n', 's') };
1053
1054         tword = getwrd(TRUE, &eol, object_codes, 1);
1055         if (tword != object_codes[0])
1056                 creaer(object_codes[0]);
1057
1058         fprintf(output, "\n         description of objects");
1059         fprintf(output, "\n                                                                density at levels");
1060         fprintf(output, "\n numb type  x-coord  y-coord x-length y-length    angle  av dens");
1061         for (i = 0; i < Spctrm.nergy; i++)
1062                 fprintf(output, " %9i", Spctrm.energy[i]);
1063
1064         // READ IN OBJECT PARAMETERS IN TEMPORARY LOCATIONS SO AS NOT TO
1065         // CLOBBER WORK AREA WHILE DOING INPUT FOR RDPICT AND RDPROJ
1066
1067         if (invoke == crea)
1068         {
1069                 f_object = new SNARK_Object_Node;
1070                 t_object = f_object;
1071         };
1072
1073         for (i = 0;; i++)
1074         {
1075                 getnxt(TRUE);           // read line of input
1076                 itype = getwrd(FALSE, &eol, &(object_codes[1]), 6);
1077
1078                 if (eol)
1079                 {
1080                         fprintf(output, "\n **** keyword elip, rect, tria, segm, sect or last is missing");
1081                         fprintf(output, "\n **** program aborted\n");
1082                         exit(-1);
1083                 };
1084
1085                 // CHECK TO SEE IF  LAST  OBJECT HAS BEEN READ
1086
1087                 if (itype == object_codes[6])
1088                         break;
1089
1090                 cxt = getnum(FALSE, &eol);
1091                 if (eol)
1092                 {
1093                         fprintf(output, "\n **** parameter cx of object %s is missing", int2str(itype));
1094                         fprintf(output, "\n **** program aborted\n");
1095                         exit(-1);
1096                 };
1097
1098                 cyt = getnum(FALSE, &eol);
1099                 if (eol)
1100                 {
1101                         fprintf(output, "\n **** parameter cy of object %s is missing", int2str(itype));
1102                         fprintf(output, "\n **** program aborted\n");
1103                         exit(-1);
1104                 };
1105
1106                 ut = getnum(FALSE, &eol);
1107                 if (eol)
1108                 {
1109                         fprintf(output, "\n **** parameter u of object %s is missing", int2str(itype));
1110                         fprintf(output, "\n **** program aborted\n");
1111                         exit(-1);
1112                 };
1113
1114                 vt = getnum(FALSE, &eol);
1115                 if (eol)
1116                 {
1117                         fprintf(output, "\n **** parameter v of object %s is missing", int2str(itype));
1118                         fprintf(output, "\n **** program aborted\n");
1119                         exit(-1);
1120                 };
1121
1122                 angt = getnum(FALSE, &eol);
1123                 if (eol)
1124                 {
1125                         fprintf(output, "\n **** parameter ang of object %s is missing", int2str(itype));
1126                         fprintf(output, "\n **** program aborted\n");
1127                         exit(-1);
1128                 };
1129
1130                 dent[0] = getnum(FALSE, &eol);
1131                 if (eol)
1132                 {
1133                         fprintf(output, "\n **** parameter den(0) of object %s is missing", int2str(itype));
1134                         fprintf(output, "\n **** program aborted\n");
1135                         exit(-1);
1136                 };
1137
1138                 // IF POLYCHROMATIC READ IN ABSORPTION FOR OTHER LEVELS
1139
1140                 if (Spctrm.nergy > 1)
1141                 {
1142                         getnxt(TRUE);           // read line of input
1143                         getwrd(FALSE, &eol, &(object_codes[7]), 1);
1144
1145                         if (eol)
1146                         {
1147                                 fprintf(output, "\n **** keyword dens is missing");
1148                                 fprintf(output, "\n **** program aborted\n");
1149                                 exit(-1);
1150                         };
1151
1152                         for (j = 1; j < Spctrm.nergy; j++)
1153                         {
1154                                 dent[j] = getnum(FALSE, &eol);
1155                                 if (eol)
1156                                 {
1157                                         fprintf(output, "\n **** parameter den(%d) of object %s is missing", j, int2str(itype));
1158                                         fprintf(output, "\n **** program aborted\n");
1159                                         exit(-1);
1160                                 }
1161                         }
1162                 }
1163
1164                 // ALSO COMPUTE THE WEIGHTED AVERAGE ABSORPTION
1165
1166                 densit = 0.0;
1167                 for (j = 0; j < Spctrm.nergy; j++)
1168                         densit += Spctrm.engwt[j] * dent[j];
1169
1170                 fprintf(output, "\n %4i %s %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f", i + 1, int2str(itype), cxt, cyt, ut, vt, angt, densit);
1171                 for (j = 0; j < Spctrm.nergy; j++)
1172                         fprintf(output, " %9.4f", dent[j]);
1173
1174                 // IDENTIFY OBJECT TYPE
1175
1176                 for (otype = SOT_elip; otype <= SOT_sect; otype++)
1177                 {
1178                         if (itype == object_codes[otype + 1])
1179                                 break;
1180                 }
1181
1182                 // IF INVOKED BY RDPICT OR RDPROJ NO NEED TO REMEMBER THE ABOVE VALUE
1183                 if (invoke == crea)
1184                 {
1185                         if ((ut < Consts.zero) || (vt < Consts.zero))
1186                                 creaer(6);
1187
1188                         // TRANSFER VALUES TO A TEMPORARY LIST
1189
1190                         t_object->data.type = (SNARK_Object_type) otype;
1191                         for (j = 0; j < Spctrm.nergy; j++)
1192                                 t_object->data.den1[j] = dent[j];
1193                         t_object->data.denav = densit;
1194                         t_object->data.cx = cxt;
1195                         t_object->data.cy = cyt;
1196                         t_object->data.u = ut;
1197                         t_object->data.v = vt;
1198                         t_object->data.ang = angt;
1199                         phi = (REAL) (angt * Consts.pi / 180.0);
1200                         t_object->data.sinang = (REAL) sin(phi);
1201                         t_object->data.cosang = (REAL) cos(phi);
1202
1203                         t_object->next = new SNARK_Object_Node;
1204                         t_object = t_object->next;
1205                 }
1206         };
1207
1208         if (invoke == crea)
1209         {
1210                 MAX_NUMBER_OF_OBJECTS = i;
1211                 objects = new SNARK_Object[MAX_NUMBER_OF_OBJECTS];
1212
1213                 t_object = f_object;
1214                 for (i = 0; i < MAX_NUMBER_OF_OBJECTS; i++)
1215                 {
1216                         objects[i].type = t_object->data.type;
1217                         for (j = 0; j < Spctrm.nergy; j++)
1218                                 objects[i].den1[j] = t_object->data.den1[j];
1219                         objects[i].denav = t_object->data.denav;
1220                         objects[i].cx = t_object->data.cx;
1221                         objects[i].cy = t_object->data.cy;
1222                         objects[i].u = t_object->data.u;
1223                         objects[i].v = t_object->data.v;
1224                         objects[i].ang = t_object->data.ang;
1225                         objects[i].sinang = t_object->data.sinang;
1226                         objects[i].cosang = t_object->data.cosang;
1227
1228                         t_object = t_object->next;
1229                         delete f_object;
1230                         f_object = t_object;
1231                 };
1232         };
1233
1234         // GRAB THE SCALE FACTOR FROM THE LAST CARD
1235
1236         scale = getnum(FALSE, &eol);
1237         if (eol)
1238         {
1239                 fprintf(output, "\n **** scale factor is missing");
1240                 fprintf(output, "\n **** program aborted\n");
1241                 exit(-1);
1242         };
1243
1244         fprintf(output,
1245                         "\n         scale factor multiplying object densities %10.4f",
1246                         scale);
1247         mean = scale;
1248
1249         if (fabs(scale) < Consts.zero)
1250                 creaer(7);
1251
1252         // IF   PICT   OR   PROJ   NO NEED TO SCALE
1253
1254         if (invoke == crea)
1255         {
1256                 // GET THE NUMBER OF OBJECTS READ IN
1257
1258                 Creacm.nobj = i;
1259                 if (Creacm.nobj != 0)
1260                 {
1261                         // SCALE DENSITIES FIRST
1262                         for (n = 0; n < Creacm.nobj; n++)
1263                         {
1264                                 objects[n].denav *= scale;
1265                                 for (j = 0; j < Spctrm.nergy; j++)
1266                                         objects[n].den1[j] *= scale;
1267                         }
1268                 }
1269         }
1270         fprintf(output, "\n");
1271
1272         //jk 11/18/2007 introducing variability
1273         //GRAB THE VARIABILITY VALUE FROM THE LAST CARD
1274         //jklukowska 12/9/8, bug 274, varseed should be read as an integer, not real number
1275         Creacm.varseed = getint(FALSE, &eol);
1276         if (eol)
1277         {
1278                 Creacm.varseed = 0;
1279         }
1280
1281         fprintf(output, "\n         seed set to %d", Creacm.varseed);
1282
1283
1284         variability = getnum(FALSE, &eol);
1285         if (eol)
1286         {
1287                 variability = 0;
1288         }
1289         else if (variability < 0)
1290         {
1291                 fprintf(output, "\n **** inhomogeneity should be >= 0");
1292                 fprintf(output, "\n **** program aborted\n");
1293                 exit(-1);
1294         }
1295
1296         if (variability != 0)
1297                 variab_set = TRUE;
1298         else if (variability == 0)
1299                 variab_set = FALSE;
1300
1301         fprintf(output, "\n         inhomogeneity set to %10.4f", variability);
1302
1303 }
1304
1305 // added by Lajos, Jan 25, 2005
1306 /*
1307  ECHOLINE PRINTS THE CONTENTS OF THE CURRENT INPUT LINE (I.E. THE ARRAY 'DATA') TO THE OUTPUT.
1308  IF flag == true, THE OUTPUT WILL START WITH '<*>'. OTHERWISE, '<#>' WILL BE USED.
1309  */
1310 void InputFile_class::echoline(BOOLEAN flag)
1311 {
1312         if (flag)
1313                 fprintf(output, "\n     <*> %s", data);
1314         else
1315                 fprintf(output, "\n     <#> %s", data);
1316 }