RH/initscatter.c at master · han-uitenbroek/RH · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
/* ------- file: -------------------------- initscatter.c -----------

       Version:       rh2.0
       Author:        Han Uitenbroek (huitenbroek@nso.edu)
       Last modified: Tue Feb 20 13:35:29 2024 --

       --------------------------                      ----------RH-- */

#include "rh.h"
#include "atom.h"
#include "atmos.h"
#include "spectrum.h"
#include "accelerate.h"
#include "error.h"
#include "statistics.h"
#include "inputs.h"


/* --- Function prototypes --                          -------------- */


/* --- Global variables --                             -------------- */

extern Atmosphere atmos;
extern Spectrum spectrum;
extern InputData input;
extern char messageStr[];


/* ------- begin -------------------------- initScatter.c ----------- */

void initScatter()
{
  const char routineName[] = "initScatter";
  register int n, kr, niter;

  double dJmax;
  Atom *atom;
  AtomicLine *line;

  /* --- Fill the radiative rates from previous solution if
         PRD lines are present --                      -------------- */

  if (atmos.NPRDactive > 0  && input.startJ == OLD_J) {
    for (n = 0;  n < atmos.Natom;  n++) {
      atom = &atmos.atoms[n];
      if (atom->active) {
	readRadRate(atom);
	for (kr = 0;  kr < atom->Nline;  kr++) {
	  line = atom->line + kr;
	  if (line->PRD) {
	    if (input.PRD_angle_dep == PRD_ANGLE_DEP)
	      PRDAngleScatter(line, LINEAR);
            else
	      PRDScatter(line, LINEAR);
	  }
	}
      }
    }
  }
  /* --- Iterate the scattering in the background. Only iterate
         when we are also planning to do main iterations -- --------- */

  if (spectrum.updateJ) {
    sprintf(messageStr,
	    " %s: Lambda iterating background scattering...\n\n",
	    routineName);
    Error(MESSAGE, routineName, messageStr);

    niter = 0;
    while (input.NmaxIter  &&  niter < input.NmaxScatter) {
      dJmax = solveSpectrum(FALSE, FALSE);
      if (dJmax < input.iterLimit) break;
      niter++;
    }
  }
  getCPU(2, TIME_POLL, "Initial solution");
}
/* ------- end ---------------------------- initScatter.c ----------- */