thermometer.cpp
1 /***************************************************************************
2  Copyright (C) 2002-2015 Kentaro Kitagawa
3  kitagawa@phys.s.u-tokyo.ac.jp
4 
5  This program is free software; you can redistribute it and/or
6  modify it under the terms of the GNU Library General Public
7  License as published by the Free Software Foundation; either
8  version 2 of the License, or (at your option) any later version.
9 
10  You should have received a copy of the GNU Library General
11  Public License and a list of authors along with this program;
12  see the files COPYING and AUTHORS.
13  ***************************************************************************/
14 //---------------------------------------------------------------------------
15 #include <math.h>
16 
17 #include "thermometer.h"
18 #include "rand.h"
19 
20 //---------------------------------------------------------------------------
21 XThermometerList::XThermometerList(const char *name, bool runtime) :
22  XCustomTypeListNode<XThermometer> (name, runtime) {
23 }
24 DECLARE_TYPE_HOLDER(XThermometerList)
25 REGISTER_TYPE(XThermometerList, LakeShore, "LakeShore")
26 ;
27 REGISTER_TYPE(XThermometerList, ScientificInstruments, "Scientific Instruments")
28 ;
29 REGISTER_TYPE(XThermometerList, ApproxThermometer, "Cubic-spline")
30 ;
31 
32 XThermometer::XThermometer(const char *name, bool runtime) :
33  XNode(name, runtime), m_tempMin(create<XDoubleNode> ("TMin", false)),
34  m_tempMax(create<XDoubleNode> ("TMax", false)) {
35  trans( *tempMin()) = 1e-3;
36  trans( *tempMax()) = 1e3;
37 }
38 
39 XLakeShore::XLakeShore(const char *name, bool runtime) :
40  XThermometer(name, runtime), m_resMin(create<XDoubleNode> ("RMin", false)),
41  m_resMax(create<XDoubleNode> ("RMax", false)), m_zu(create<
42  XDoubleListNode> ("Zu", false)), m_zl(create<XDoubleListNode> (
43  "Zl", false)), m_ai(create<XDouble2DNode> ("Ai", false)) {
44 }
45 
46 double XLakeShore::getRawValue(double temp) const {
47  Snapshot shot( *this);
48  //using Newton's method
49  double x, y, dypdx, val;
50  if(temp < shot[ *tempMin()])
51  return shot[ *resMax()];
52  if(temp > shot[ *tempMax()])
53  return shot[ *resMin()];
54  // x = (log10(RMin) + log10(RMax)) / 2;
55  val = shot[ *resMin()];
56  for(double dy = 0.0001;; dy *= 2) {
57  if (dy > 1.0)
58  return shot[ *resMin()];
59  double t = randMT19937();
60  x = (log10(shot[ *resMax()]) * t + log10(shot[ *resMin()]) * (1 - t));
61  for(int i = 0; i < 100; i++) {
62  y = getTemp(pow(10, x)) - temp;
63  if(fabs(y) < dy) {
64  val = pow(10, x);
65  return val;
66  }
67  dypdx = (y - (getTemp(pow(10, x - 0.00001)) - temp)) / 0.00001;
68  if(dypdx != 0)
69  x -= y / dypdx;
70  if((x > log10(shot[ *resMax()])) || (x < log10(shot[ *resMin()]))
71  || (dypdx == 0)) {
72  double t = randMT19937();
73  x = (log10(shot[ *resMax()]) * t + log10(shot[ *resMin()]) * (1 - t));
74  }
75  }
76  }
77  return val;
78 }
79 
80 double XLakeShore::getTemp(double res) const {
81  Snapshot shot( *this);
82  double temp = 0, z, u = 0;
83  if(res > shot[ *resMax()])
84  return shot[ *tempMin()];
85  if(res < shot[ *resMin()])
86  return shot[ *tempMax()];
87  z = log10(res);
88  unsigned int n;
89  if( !shot.size(zu()))
90  return 0;
91  const auto &zu_list( *shot.list(zu()));
92  if( !shot.size(zl()))
93  return 0;
94  const auto &zl_list( *shot.list(zl()));
95  for(n = 0; n < zu_list.size(); n++) {
96  double zu = shot[ *static_pointer_cast<XDoubleNode> (zu_list.at(n))];
97  double zl = shot[ *static_pointer_cast<XDoubleNode> (zl_list.at(n))];
98  u = (z - zu + z - zl) / (zu - zl);
99  if ((u >= -1) && (u <= 1))
100  break;
101  }
102  if(n >= zu_list.size())
103  return 0;
104  if( !shot.size(ai()))
105  return 0;
106  const auto &ai_list( *shot.list(ai()));
107  if( !shot.size(ai_list[n]))
108  return 0;
109  const auto &ai_n_list( *shot.list(ai_list[n]));
110  for (unsigned int i = 0; i < ai_n_list.size(); i++) {
111  double ai_n_i = shot[ *static_pointer_cast<XDoubleNode> (ai_n_list.at(i))];
112  temp += ai_n_i * cos(i * acos(u));
113  }
114  return temp;
115 }
116 
117 XScientificInstruments::XScientificInstruments(const char *name, bool runtime) :
118  XThermometer(name, runtime), m_resMin(create<XDoubleNode> ("RMin", false)),
119  m_resMax(create<XDoubleNode> ("RMax", false)), m_abcde(create<
120  XDoubleListNode> ("ABCDE", false)), m_abc(create<XDoubleListNode> (
121  "ABC", false)), m_rCrossover(create<XDoubleNode> ("RCrossover",
122  false)) {
123 }
124 
125 double XScientificInstruments::getRawValue(double temp) const {
126  Snapshot shot( *this);
127  //using Newton's method
128  double x, y, dypdx, val;
129  if(temp < shot[ *tempMin()])
130  return shot[ *resMax()];
131  if(temp > shot[ *tempMax()])
132  return shot[ *resMin()];
133  // x = (log10(RMin) + log10(RMax)) / 2;
134  val = shot[ *resMin()];
135  for(double dy = 0.0001;; dy *= 2) {
136  if(dy > 1.0)
137  return shot[ *resMin()];
138  double t = randMT19937();
139  x = (log10(shot[ *resMax()]) * t + log10(shot[ *resMin()]) * (1 - t));
140  for(int i = 0; i < 100; i++) {
141  y = getTemp(pow(10, x)) - temp;
142  if(fabs(y) < dy) {
143  val = pow(10, x);
144  return val;
145  }
146  dypdx = (y - (getTemp(pow(10, x - 0.00001)) - temp)) / 0.00001;
147  if(dypdx != 0)
148  x -= y / dypdx;
149  if((x > log10(shot[ *resMax()])) || (x < log10(shot[ *resMin()]))
150  || (dypdx == 0)) {
151  double t = randMT19937();
152  x = (log10(shot[ *resMax()]) * t + log10(shot[ *resMin()]) * (1 - t));
153  }
154  }
155  }
156  return val;
157 }
158 
159 double XScientificInstruments::getTemp(double res) const {
160  Snapshot shot( *this);
161  if(res > shot[ *resMax()])
162  return shot[ *tempMin()];
163  if(res < shot[ *resMin()])
164  return shot[ *tempMax()];
165  double y = 0.0;
166  double lx = log(res);
167  if(res > shot[ *rCrossover()]) {
168  if( !shot.size(abcde())) return 0;
169  const auto &abcde_list( *shot.list(abcde()));
170  if(abcde_list.size() >= 5) {
171  double a = shot[ *static_pointer_cast<XDoubleNode> (abcde_list.at(0))];
172  double b = shot[ *static_pointer_cast<XDoubleNode> (abcde_list.at(1))];
173  double c = shot[ *static_pointer_cast<XDoubleNode> (abcde_list.at(2))];
174  double d = shot[ *static_pointer_cast<XDoubleNode> (abcde_list.at(3))];
175  double e = shot[ *static_pointer_cast<XDoubleNode> (abcde_list.at(4))];
176  y = (a + c * lx + e * lx * lx) / (1.0 + b * lx + d * lx * lx);
177  }
178  return y;
179  } else {
180  if( !shot.size(abc())) return 0;
181  const auto &abc_list( *shot.list(abc()));
182  if(abc_list.size() >= 3) {
183  double a = shot[ *static_pointer_cast<XDoubleNode> (abc_list.at(0))];
184  double b = shot[ *static_pointer_cast<XDoubleNode> (abc_list.at(1))];
185  double c = shot[ *static_pointer_cast<XDoubleNode> (abc_list.at(2))];
186  y = 1.0 / (a + b * res * lx + c * res * res);
187  }
188  return y;
189  }
190 }
191 
192 XApproxThermometer::XApproxThermometer(const char *name, bool runtime) :
193  XThermometer(name, runtime), m_resList(create<XDoubleListNode> ("ResList",
194  false)), m_tempList(create<XDoubleListNode> ("TempList", false)) {
195 }
196 
197 double XApproxThermometer::getTemp(double res) const {
198  local_shared_ptr<CSplineApprox> approx(m_approx);
199  Snapshot shot( *this);
200  if( !approx) {
201  std::map<double, double> pts;
202  if( !shot.size(m_resList))
203  return 0;
204  const XNode::NodeList &res_list( *shot.list(m_resList));
205  if( !shot.size(m_tempList))
206  return 0;
207  const auto &temp_list( *shot.list(m_tempList));
208  for(unsigned int i = 0; i < std::min(res_list.size(), temp_list.size()); i++) {
209  double res = shot[ *static_pointer_cast<XDoubleNode> (res_list.at(i))];
210  double temp = shot[ *static_pointer_cast<XDoubleNode> (temp_list.at(i))];
211  pts.insert(std::pair<double, double>(log(res), log(temp)));
212  }
213  if(pts.size() < 4)
214  throw XKameError(i18n(
215  "XApproxThermometer, Too small number of points"), __FILE__, __LINE__);
216  approx.reset(new CSplineApprox(pts));
217  m_approx = approx;
218  }
219  return exp(approx->approx(log(res)));
220 }
221 double XApproxThermometer::getRawValue(double temp) const {
222  local_shared_ptr<CSplineApprox> approx(m_approx_inv);
223  Snapshot shot( *this);
224  if( !approx) {
225  std::map<double, double> pts;
226  if( !shot.size(m_resList))
227  return 0;
228  const auto &res_list( *shot.list(m_resList));
229  if( !shot.size(m_tempList))
230  return 0;
231  const auto &temp_list( *shot.list(m_tempList));
232  for(unsigned int i = 0; i < std::min(res_list.size(), temp_list.size()); i++) {
233  double res = shot[ *static_pointer_cast<XDoubleNode> (res_list.at(i))];
234  double temp = shot[ *static_pointer_cast<XDoubleNode> (temp_list.at(i))];
235  pts.insert(std::pair<double, double>(log(temp), log(res)));
236  }
237  if(pts.size() < 4)
238  throw XKameError(i18n(
239  "XApproxThermometer, Too small number of points"), __FILE__, __LINE__);
240  approx.reset(new CSplineApprox(pts));
241  m_approx_inv = approx;
242  }
243  return exp(approx->approx(log(temp)));
244 }
245 
Generated for KAME4 by  doxygen 1.8.3