Line data Source code
1 : /**
2 : Bojan Nikolic <b.nikolic@mrao.cam.ac.uk>, <bojan@bnikolic.co.uk>
3 : Initial version January 2010.
4 : Maintained by ESO since 2013. Incorporating code originally in
5 : apps/almaabs_i.cpp.
6 :
7 : This file is part of LibAIR and is licensed under GNU Public
8 : License Version 2
9 :
10 : \file dtdltools.cpp
11 : Renamed dtdltools.cc 2023
12 : */
13 :
14 : #include <iostream>
15 : #include <cmath>
16 : #include <array>
17 :
18 : #include "bnmin1/src/nestedsampler.h"
19 :
20 : #include "dtdltools.h"
21 : #include "model_make.h"
22 : #include "numalgo.h"
23 :
24 :
25 : namespace LibAIR2 {
26 :
27 0 : void dTdLMom1(const std::list<Minim::WPPoint> &l,
28 : Minim::ModelDesc &md,
29 : const WVRAtmoQuants &model,
30 : double Z,
31 : double thresh,
32 : double *res)
33 : {
34 0 : std::vector<double> scratch(4, 0.0);
35 0 : const double Pthresh= Z*thresh;
36 0 : for(size_t i=0; i<4; ++i)
37 0 : res[i]=0;
38 :
39 0 : for(std::list<Minim::WPPoint>::const_iterator i=l.begin();
40 0 : i!= l.end();
41 0 : ++i)
42 : {
43 0 : const double w=i->w *exp(- i->ll);
44 0 : if (w > Pthresh)
45 : {
46 0 : md.put(i->p);
47 0 : model.dTdL_ND(scratch);
48 0 : for(size_t j=0; j<scratch.size(); ++j)
49 0 : res[j]+=(scratch[j]*w);
50 : }
51 : }
52 0 : if(Z==0.)
53 : {
54 0 : std::cout << "Error: Cannot calculate dTdL Moment 1, evidence is zero." << std::endl;
55 0 : std::cerr << "Error: Cannot calculate dTdL Moment 1, evidence is zero." << std::endl;
56 : }
57 : else
58 : {
59 0 : for(size_t j=0; j<scratch.size(); ++j)
60 : {
61 0 : res[j]/=Z;
62 : }
63 : }
64 0 : }
65 :
66 0 : void dTdLMom2(const std::list<Minim::WPPoint> &l,
67 : Minim::ModelDesc &md,
68 : const WVRAtmoQuants &model,
69 : const double *m1,
70 : double Z,
71 : double thresh,
72 : double *res
73 : )
74 : {
75 0 : std::vector<double> scratch(4, 0.0);
76 0 : const double Pthresh= Z*thresh;
77 0 : for(size_t i=0; i<4; ++i)
78 0 : res[i]=0;
79 :
80 0 : for(std::list<Minim::WPPoint>::const_iterator i=l.begin();
81 0 : i!= l.end();
82 0 : ++i)
83 : {
84 0 : const double w=i->w *exp(- i->ll);
85 0 : if (w > Pthresh)
86 : {
87 0 : md.put(i->p);
88 0 : model.dTdL_ND(scratch);
89 0 : for(size_t j=0; j<scratch.size(); ++j)
90 : {
91 0 : res[j]+=std::pow(scratch[j]-m1[j],2)*w;
92 : }
93 : }
94 : }
95 0 : if(Z==0.)
96 : {
97 0 : std::cout << "Error: Cannot calculate dTdL Moment 2, evidence is zero." << std::endl;
98 0 : std::cerr << "Error: Cannot calculate dTdL Moment 2, evidence is zero." << std::endl;
99 : }
100 : else
101 : {
102 0 : for(size_t j=0; j<scratch.size(); ++j)
103 : {
104 0 : res[j]/=Z;
105 : }
106 : }
107 0 : }
108 :
109 : struct CenFD_bind {
110 : WVRAtmoQuantModel &md;
111 0 : CenFD_bind(WVRAtmoQuantModel &md_) : md(md_) { }
112 0 : std::array<double, 4> operator( ) (double d) { return md.evalFn(d,"n"); }
113 : };
114 :
115 0 : void dTdL2_ND(WVRAtmoQuantModel &m,
116 : std::vector<double> &res)
117 : {
118 0 : Minim::ModelDesc md(m);
119 0 : CenFD_bind bind_instance(m);
120 : std::array< std::array<double, 4> ,3> r=
121 0 : CenFDV<3,4>(bind_instance,
122 0 : *md["n"]->p,
123 : 0.001);
124 0 : const double cv=std::pow(SW_WaterToPath_Simplified(1.0,
125 0 : *md["T"]->p),
126 : -2);
127 0 : res=std::vector<double>(r[2].begin(), r[2].end());
128 0 : for(double &x: res)
129 : {
130 0 : x*=cv;
131 : }
132 :
133 0 : }
134 :
135 0 : void dTdL3_ND(WVRAtmoQuantModel &m,
136 : std::vector<double> &res)
137 : {
138 0 : Minim::ModelDesc md(m);
139 0 : CenFD_bind bind_instance(m);
140 : std::array< std::array<double, 4> ,4> r=
141 0 : CenFDV<4,4>(bind_instance,
142 0 : *md["n"]->p,
143 : 0.001);
144 0 : const double cv=std::pow(SW_WaterToPath_Simplified(1.0,
145 0 : *md["T"]->p),
146 : -3);
147 0 : res=std::vector<double>(r[3].begin(), r[3].end());
148 0 : for(double &x: res)
149 : {
150 0 : x*=cv;
151 : }
152 :
153 0 : }
154 :
155 0 : void dTdTAtm(WVRAtmoQuantModel &m,
156 : std::vector<double> &res)
157 : {
158 0 : Minim::ModelDesc md(m);
159 0 : CenFD_bind bind_instance(m);
160 : std::array< std::array<double, 4> ,2> r=
161 0 : CenFDV<2,4>(bind_instance,
162 0 : *md["T"]->p,
163 : 0.01);
164 0 : res=std::vector<double>(r[1].begin(), r[1].end());
165 0 : }
166 :
167 : }
168 :
169 :
|
