-
Notifications
You must be signed in to change notification settings - Fork 17
/
redistribute.c
117 lines (90 loc) · 3.21 KB
/
redistribute.c
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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
/* ------- file: -------------------------- redistribute.c ----------
Version: rh2.0
Author: Han Uitenbroek ([email protected])
Last modified: Wed Apr 1 14:01:22 2009 --
-------------------------- ----------RH-- */
/* --- Administers iterations to redistribute intensity in PRD line while
keeping the population number fixed. -- -------------- */
#include <stdlib.h>
#include <math.h>
#include "rh.h"
#include "atom.h"
#include "atmos.h"
#include "spectrum.h"
#include "accelerate.h"
#include "error.h"
#include "inputs.h"
/* --- Function prototypes -- -------------- */
/* --- Global variables -- -------------- */
extern Atmosphere atmos;
extern Spectrum spectrum;
extern InputData input;
extern CommandLine commandline;
extern char messageStr[];
/* ------- begin -------------------------- Redistribute.c ---------- */
void Redistribute(int NmaxIter, double iterLimit)
{
const char routineName[] = "Redistribute";
register int kr, nact;
bool_t quiet, accel, eval_operator, redistribute;
enum Interpolation representation;
int niter, Nlamu;
double drho, drhomax, drhomaxa;
Atom *atom;
AtomicLine *line;
for (nact = 0; nact < atmos.Nactiveatom; nact++) {
atom = atmos.activeatoms[nact];
/* --- Initialize structures for Ng acceleration PRD iteration -- */
for (kr = 0; kr < atom->Nline; kr++) {
line = &atom->line[kr];
if (line->PRD && line->Ng_prd == NULL) {
if (input.PRD_angle_dep == PRD_ANGLE_DEP)
Nlamu = 2*atmos.Nrays * line->Nlambda * atmos.Nspace;
else
Nlamu = line->Nlambda*atmos.Nspace;
line->Ng_prd = NgInit(Nlamu, input.PRD_Ngdelay,
input.PRD_Ngorder, input.PRD_Ngperiod,
line->rho_prd[0]);
}
}
}
/* --- Iterate over scattering integral while keeping populations
fixed -- -------------- */
niter = 1;
while (niter <= NmaxIter) {
drhomaxa = 0.0;
for (nact = 0; nact < atmos.Nactiveatom; nact++) {
atom = atmos.activeatoms[nact];
drhomax = 0.0;
for (kr = 0; kr < atom->Nline; kr++) {
line = &atom->line[kr];
if (line->PRD) {
switch (input.PRD_angle_dep) {
case PRD_ANGLE_INDEP:
PRDScatter(line, representation=LINEAR);
break;
case PRD_ANGLE_APPROX:
PRDAngleApproxScatter(line, representation=LINEAR);
break;
case PRD_ANGLE_DEP:
PRDAngleScatter(line, representation=LINEAR);
break;
}
accel = Accelerate(line->Ng_prd, line->rho_prd[0]);
sprintf(messageStr, " PRD: iter #%d, atom %s, line %d,",
line->Ng_prd->count-1, atom->ID, kr);
drho = MaxChange(line->Ng_prd, messageStr, quiet=FALSE);
sprintf(messageStr, (accel) ? " (accelerated)\n" : "\n");
Error(MESSAGE, routineName, messageStr);
drhomax = MAX(drho, drhomax);
}
drhomaxa = MAX(drhomax, drhomaxa);
}
}
/* --- Solve transfer equation with fixed populations -- -------- */
solveSpectrum(eval_operator=FALSE, redistribute=TRUE);
if (drhomaxa < iterLimit) break;
niter++;
}
}
/* ------- end ---------------------------- Redistribute.c ---------- */