v0.14.0
Loading...
Searching...
No Matches
CrudeClimateModel.hpp
Go to the documentation of this file.
1/*
2 * This file is part of MoFEM.
3 * MoFEM is free software: you can redistribute it and/or modify it under
4 * the terms of the GNU Lesser General Public License as published by the
5 * Free Software Foundation, either version 3 of the License, or (at your
6 * option) any later version.
7 *
8 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
9 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
11 * License for more details.
12 *
13 * You should have received a copy of the GNU Lesser General Public
14 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
15
16#ifndef __GROUND_TIME_DATA_HPP
17#define __GROUND_TIME_DATA_HPP
18
19//#include <boost/program_options.hpp>
20//#include <iostream>
21//#include <fstream>
22//#include <iterator>
23
24using namespace std;
25namespace po = boost::program_options;
26
27/**
28 *
29 * This is based on:
30 * NEW MODELS FOR SIMULATING DAILY MINIMUM, DAILY MAXIMUM AND HOURLY OUTDOOR TEMPERATURES
31 * H Bulut1, O. Buyukalaca2 and T. Yılmaz2
32 *
33 * However some modifications are introduced, that model can be set using data
34 * easy found on internet.
35 *
36 */
38
39 CrudeClimateModel(const char file[]): GenericClimateModel() {
40
41 ierr = setDefault(); CHKERRABORT(PETSC_COMM_WORLD,ierr);
42 ierr = readFile(file); CHKERRABORT(PETSC_COMM_WORLD,ierr);
43 }
44
47 int Low,High;
48
49 double fractinalDayTime(double t) {
50 const double nb_days = 365;
51 double n = t/(60*60*24);
52 return n - floor(n/nb_days)*nb_days;
53 }
54 double wetherTime(double n) {
55 const double nb_days = 365;
56 double a = -0.1e1 / (double) nb_days / (double) (High - Low - nb_days) * M_PI * (double) (2 * High - 2 * Low - nb_days) / (double) (High - Low);
57 double b = M_PI * (double) (2 * High * High - 2 * Low * Low - 2 * nb_days * Low - nb_days * nb_days) / (double) nb_days / (double) (High - Low - nb_days) / (double) (High - Low);
58 double c = -0.1e1 / (double) nb_days / (double) (High - Low - nb_days) * M_PI * (double) (2 * High - 2 * Low - nb_days) / (double) (High - Low);
59 return a+b*n+c*n*n;
60 }
61 double funTday(double t) {
62 double n = fractinalDayTime(t);
64 }
65 double funTnigth(double t) {
66 double n = fractinalDayTime(t);
68 }
69
70 double Td; // dwe point
71 struct tm startTime;
72 time_t t0;
73
74
75 PetscErrorCode readFile(const char file[]) {
76 PetscFunctionBegin;
77 try {
78 ifstream ini_file(file);
79 po::variables_map vm;
80 po::options_description config_file_options;
81 config_file_options.add_options()
82 ("longitude",po::value<double>(&spaData.longitude)->default_value(0.127))
83 ("latitude",po::value<double>(&spaData.latitude)->default_value(51.5072))
84 ("year",po::value<int>(&spaData.year)->default_value(2014))
85 ("month",po::value<int>(&spaData.month)->default_value(6))
86 ("day",po::value<int>(&spaData.day)->default_value(21))
87 ("hour",po::value<int>(&spaData.hour)->default_value(0))
88 ("minute",po::value<int>(&spaData.minute)->default_value(0))
89 ("second",po::value<double>(&spaData.second)->default_value(0))
90 ("timezone",po::value<double>(&spaData.timezone)->default_value(0))
91 ("elevation",po::value<double>(&spaData.elevation)->default_value(10))
92 ("TdayAtSummer",po::value<double>(&TdaySummer)->default_value(10),"Day temerature at summer")
93 ("TnigthAtSummer",po::value<double>(&TnigthSummer)->default_value(8),"Nigth temerature at summer")
94 ("TdayAtWinter",po::value<double>(&TdayWinter)->default_value(-8),"Day temerature at winter")
95 ("TnigthAtWinter",po::value<double>(&TnigthWinter)->default_value(-10),"Nigth temerature at winter")
96 ("DayOfLowTemperature",po::value<int>(&Low)->default_value(0))
97 ("DayOfHighTemperature",po::value<int>(&High)->default_value(182))
98 ("DewPoint",po::value<double>(&Td)->default_value(0),"Dew point in Celsius degrees")
99 ("u10",po::value<double>(&u10)->default_value(2.7))
100 ("CR",po::value<double>(&CR)->default_value(0))
101 ("Pressure",po::value<double>(&P)->default_value(101325))
102 ("Rs",po::value<double>(&Rs)->default_value(0));
103 store(parse_config_file(ini_file,config_file_options), vm);
104 po::notify(vm);
105
106 startTime.tm_sec = spaData.second;
107 startTime.tm_min = spaData.minute;
108 startTime.tm_hour = spaData.hour;
109 startTime.tm_mday = spaData.day;
110 startTime.tm_mon = spaData.month-1;
111 startTime.tm_year = spaData.year-1900;
112 startTime.tm_isdst = 0;
113 t0 = mktime(&startTime);
114
115 } catch (const std::exception& ex) {
116 ostringstream ss;
117 ss << ex.what() << endl;
118 SETERRQ(PETSC_COMM_SELF,MOFEM_STD_EXCEPTION_THROW,ss.str().c_str());
119 }
120
121 PetscFunctionReturn(0);
122 }
123
124 double funT(double t,double t0 = 0) {
125 double Tmax = funTday(t);
126 double Tmin = funTnigth(t);
127
128 //cerr << "Tmax " << Tmax << " Tmin " << Tmin << endl;
129
130 //spaData.suntransit,spaData.sunrise,spaData.sunset
131 double td = (spaData.sunset-spaData.sunrise);
132
133 double tmin = 12-td/2;
134 double tmax = 12+(tmin*(12-tmin))/13.5;
135
136 double h = (t-floor(t/(3600*24))*3600*24)/3600.;
137
138 double T;
139 if((h<tmin)||(h>tmax)) {
140 double h_ = h;
141 if(h<tmin) h_+=24;
142 T = (Tmax-(Tmax-Tmin)*pow(sin(0.5*M_PI*(h_-tmax)/(24+tmin-tmax)),1.2));
143 } else {
144 T = Tmin+(Tmax-Tmin)*pow(sin(0.5*M_PI*(h-tmin)/(tmax-tmin)),1.4);
145 }
146
147 return T;
148 }
149
150 PetscErrorCode testCode(double duration,double step) {
151 PetscFunctionBegin;
152
153 /*cerr << spaData.second << endl;
154 cerr << spaData.minute << endl;
155 cerr << spaData.hour << endl;
156 cerr << spaData.day << endl;
157 cerr << spaData.month-1 << endl;
158 cerr << spaData.year << endl;*/
159 //cerr << TdaySummer << endl;
160 //cerr << TdayWinter << endl;
161 //cerr << TnigthSummer << endl;
162 //cerr << TnigthWinter << endl;
163
164 #ifdef BOOST_IOSTREAMS_TEE_HPP_INCLUDED
165 typedef tee_device<ostream, ofstream> TeeDevice;
166 typedef stream<TeeDevice> TeeStream;
167 ofstream ofs("time_variation_model.txt");
168 TeeDevice my_tee(cout, ofs);
169 TeeStream my_split(my_tee);
170 #endif //BOOST_IOSTREAMS_TEE_HPP_INCLUDED
171
172 time_t t = t0;
173 for(;t<t0+duration;t+=step) {
174
175 struct tm current_time;
176 current_time = *gmtime(&t);
177 spaData.second = current_time.tm_sec;
178 spaData.minute = current_time.tm_min;
179 spaData.hour = current_time.tm_hour;
180 spaData.day = current_time.tm_mday;
181 spaData.month = current_time.tm_mon+1;
182 spaData.year = current_time.tm_year+1900;
183
185 int r;
187 if(r) {
188 SETERRQ1(PETSC_COMM_SELF,1,"wrong input data for solar position calulator error codde %d",r);
189 }
190
191 double T = funT(t,t0);
192
193 #ifndef BOOST_IOSTREAMS_TEE_HPP_INCLUDED
194 PetscPrintf(PETSC_COMM_WORLD,"%3.4f %3.4f\n",((double)t-(double)t0)/(60*60*24),T);
195 #else //BOOST_IOSTREAMS_TEE_HPP_INCLUDED
196 my_split.precision(3);
197 my_split << ((double)t-(double)t0)/(60*60*24) << " " << T << endl;
198 #endif //BOOST_IOSTREAMS_TEE_HPP_INCLUDED
199
200 }
201
202 PetscFunctionReturn(0);
203
204 }
205
206 PetscErrorCode setDefault() {
207 PetscFunctionBegin;
208
209 spaData.delta_ut1 = 0; // Fractional second difference between UTC and UT which is used
210 // to adjust UTC for earth's irregular rotation rate and is derived
211 // from observation only and is reported in this bulletin:
212 // http://maia.usno.navy.mil/ser7/ser7.dat,
213 // where delta_ut1 = DUT1
214 // valid range: -1 to 1 second (exclusive), error code 17
215
216 spaData.delta_t = 0; // Difference between earth rotation time and terrestrial time
217 // It is derived from observation only and is reported in this
218 // bulletin: http://maia.usno.navy.mil/ser7/ser7.dat,
219 // where delta_t = 32.184 + (TAI-UTC) - DUT1
220 // valid range: -8000 to 8000 seconds, error code: 7
221
222 spaData.pressure = 1013.25; // Annual average local pressure [millibars]
223 // valid range: 0 to 5000 millibars, error code: 12
224
225 spaData.temperature = 20; // Annual average local temperature [degrees Celsius]
226 // valid range: -273 to 6000 degrees Celsius, error code; 13
227
228
229 spaData.slope = 0; // Surface slope (measured from the horizontal plane)
230 // valid range: -360 to 360 degrees, error code: 14
231
232 spaData.azm_rotation = 0; // Surface azimuth rotation (measured from south to projection of
233 // surface normal on horizontal plane, negative east)
234 // valid range: -360 to 360 degrees, error code: 15
235
236 spaData.atmos_refract = 0.5667; // Atmospheric refraction at sunrise and sunset (0.5667 deg is typical)
237 // valid range: -5 to 5 degrees, error code: 16
238
239 PetscFunctionReturn(0);
240 }
241
243
244 /** \brief set time in days
245 */
246 PetscErrorCode set(double time = 0) {
247 PetscFunctionBegin;
248
249 time_t t = t0;
250 t += time*60*60*24;
251
252 struct tm current_time;
253 current_time = *gmtime(&t);
254 spaData.second = current_time.tm_sec;
255 spaData.minute = current_time.tm_min;
256 spaData.hour = current_time.tm_hour;
257 spaData.day = current_time.tm_mday;
258 spaData.month = current_time.tm_mon+1;
259 spaData.year = current_time.tm_year+1900;
260
262 int r;
264 if(r) {
265 SETERRQ1(PETSC_COMM_SELF,1,"wrong input data for solar position calulator error codde %d",r);
266 }
269
270 Ta = funT(t,t0); // sets air temperature
271
272 PetscPrintf(PETSC_COMM_WORLD,
273 "Suntransit %3.2f Sunrise %3.2f Sunset %3.2f Ta = %3.2f \n" ,
275 Ta);
276 PetscPrintf(PETSC_COMM_WORLD,"Time %s",asctime(&current_time));
277
278
279 PetscFunctionReturn(0);
280 }
281
282};
283
284
285#endif //__GROUND_TIME_DATA_HPP
constexpr double a
static PetscErrorCode ierr
@ MOFEM_STD_EXCEPTION_THROW
Definition definitions.h:39
FTensor::Index< 'n', SPACE_DIM > n
tee_device< std::ostream, std::ofstream > TeeDevice
const double c
speed of light (cm/ns)
const double T
double h
constexpr double t
plate stiffness
Definition plate.cpp:59
int spa_calculate(spa_data *spa)
Definition spa.c:1131
@ SPA_ZA_RTS
Definition spa.h:60
double wetherTime(double n)
double funTday(double t)
PetscErrorCode readFile(const char file[])
PetscErrorCode testCode(double duration, double step)
CrudeClimateModel(const char file[])
double fractinalDayTime(double t)
PetscErrorCode setDefault()
PetscErrorCode set(double time=0)
set time in days
double funT(double t, double t0=0)
double funTnigth(double t)
Definition spa.h:65
double zenith
Definition spa.h:171
double suntransit
Definition spa.h:176
double second
Definition spa.h:73
double latitude
Definition spa.h:94
int function
Definition spa.h:116
double azimuth
Definition spa.h:173
int hour
Definition spa.h:71
double timezone
Definition spa.h:88
double longitude
Definition spa.h:91
double pressure
Definition spa.h:100
double sunrise
Definition spa.h:177
double sunset
Definition spa.h:178
int day
Definition spa.h:70
double azm_rotation
Definition spa.h:109
double elevation
Definition spa.h:97
double temperature
Definition spa.h:103
int minute
Definition spa.h:72
double delta_ut1
Definition spa.h:75
double slope
Definition spa.h:106
double delta_t
Definition spa.h:82
int month
Definition spa.h:69
double atmos_refract
Definition spa.h:113
int year
Definition spa.h:68