covid-sim
Dist.cpp
1 #include <stdlib.h>
2 #include <math.h>
3 
4 #include "Constants.h"
5 #include "Dist.h"
6 #include "Param.h"
7 
8 double sinx[361], cosx[361], asin2sqx[1001];
9 
13 
14 double dist2UTM(double x1, double y1, double x2, double y2)
15 {
16  double x, y, cy1, cy2, yt, xi, yi;
17 
18  x = fabs(x1 - x2) / 2;
19  y = fabs(y1 - y2) / 2;
20  xi = floor(x);
21  yi = floor(y);
22  x -= xi;
23  y -= yi;
24  x = (1 - x) * sinx[(int)xi] + x * sinx[((int)xi) + 1];
25  y = (1 - y) * sinx[(int)yi] + y * sinx[((int)yi) + 1];
26  yt = fabs(y1 + P.SpatialBoundingBox[1]);
27  yi = floor(yt);
28  cy1 = yt - yi;
29  cy1 = (1 - cy1) * cosx[((int)yi)] + cy1 * cosx[((int)yi) + 1];
30  yt = fabs(y2 + P.SpatialBoundingBox[1]);
31  yi = floor(yt);
32  cy2 = yt - yi;
33  cy2 = (1 - cy2) * cosx[((int)yi)] + cy2 * cosx[((int)yi) + 1];
34  x = fabs(1000 * (y * y + x * x * cy1 * cy2));
35  xi = floor(x);
36  x -= xi;
37  y = (1 - x) * asin2sqx[((int)xi)] + x * asin2sqx[((int)xi) + 1];
38  return 4 * EARTHRADIUS * EARTHRADIUS * y;
39 }
40 double dist2(Person* a, Person* b)
41 {
42  double x, y;
43 
44  if (P.DoUTM_coords)
45  return dist2UTM(Households[a->hh].loc_x, Households[a->hh].loc_y, Households[b->hh].loc_x, Households[b->hh].loc_y);
46  else
47  {
48  x = fabs(Households[a->hh].loc_x - Households[b->hh].loc_x);
49  y = fabs(Households[a->hh].loc_y - Households[b->hh].loc_y);
50  if (P.DoPeriodicBoundaries)
51  {
52  if (x > P.in_degrees_.width_ * 0.5) x = P.in_degrees_.width_ - x;
53  if (y > P.in_degrees_.height_ * 0.5) y = P.in_degrees_.height_ - y;
54  }
55  return x * x + y * y;
56  }
57 }
58 double dist2_cc(Cell* a, Cell* b)
59 {
60  double x, y;
61  int l, m;
62 
63  l = (int)(a - Cells);
64  m = (int)(b - Cells);
65  if (P.DoUTM_coords)
66  return dist2UTM(P.in_cells_.width_ * fabs((double)(l / P.nch)), P.in_cells_.height_ * fabs((double)(l % P.nch)),
67  P.in_cells_.width_ * fabs((double)(m / P.nch)), P.in_cells_.height_ * fabs((double)(m % P.nch)));
68  else
69  {
70  x = P.in_cells_.width_ * fabs((double)(l / P.nch - m / P.nch));
71  y = P.in_cells_.height_ * fabs((double)(l % P.nch - m % P.nch));
72  if (P.DoPeriodicBoundaries)
73  {
74  if (x > P.in_degrees_.width_ * 0.5) x = P.in_degrees_.width_ - x;
75  if (y > P.in_degrees_.height_ * 0.5) y = P.in_degrees_.height_ - y;
76  }
77  return x * x + y * y;
78  }
79 }
80 double dist2_cc_min(Cell* a, Cell* b)
81 {
82  double x, y;
83  int l, m, i, j;
84 
85  l = (int)(a - Cells);
86  m = (int)(b - Cells);
87  i = l; j = m;
88  if (P.DoUTM_coords)
89  {
90  if (P.in_cells_.width_ * ((double)abs(m / P.nch - l / P.nch)) > PI)
91  {
92  if (m / P.nch > l / P.nch)
93  j += P.nch;
94  else if (m / P.nch < l / P.nch)
95  i += P.nch;
96  }
97  else
98  {
99  if (m / P.nch > l / P.nch)
100  i += P.nch;
101  else if (m / P.nch < l / P.nch)
102  j += P.nch;
103  }
104  if (m % P.nch > l % P.nch)
105  i++;
106  else if (m % P.nch < l % P.nch)
107  j++;
108  return dist2UTM(P.in_cells_.width_ * fabs((double)(i / P.nch)), P.in_cells_.height_ * fabs((double)(i % P.nch)),
109  P.in_cells_.width_ * fabs((double)(j / P.nch)), P.in_cells_.height_ * fabs((double)(j % P.nch)));
110  }
111  else
112  {
113  if ((P.DoPeriodicBoundaries) && (P.in_cells_.width_ * ((double)abs(m / P.nch - l / P.nch)) > P.in_degrees_.width_ * 0.5))
114  {
115  if (m / P.nch > l / P.nch)
116  j += P.nch;
117  else if (m / P.nch < l / P.nch)
118  i += P.nch;
119  }
120  else
121  {
122  if (m / P.nch > l / P.nch)
123  i += P.nch;
124  else if (m / P.nch < l / P.nch)
125  j += P.nch;
126  }
127  if ((P.DoPeriodicBoundaries) && (P.in_degrees_.height_ * ((double)abs(m % P.nch - l % P.nch)) > P.in_degrees_.height_ * 0.5))
128  {
129  if (m % P.nch > l % P.nch)
130  j++;
131  else if (m % P.nch < l % P.nch)
132  i++;
133  }
134  else
135  {
136  if (m % P.nch > l % P.nch)
137  i++;
138  else if (m % P.nch < l % P.nch)
139  j++;
140  }
141  x = P.in_cells_.width_ * fabs((double)(i / P.nch - j / P.nch));
142  y = P.in_cells_.height_ * fabs((double)(i % P.nch - j % P.nch));
143  if (P.DoPeriodicBoundaries)
144  {
145  if (x > P.in_degrees_.width_ * 0.5) x = P.in_degrees_.width_ - x;
146  if (y > P.in_degrees_.height_ * 0.5) y = P.in_degrees_.height_ - y;
147  }
148  return x * x + y * y;
149  }
150 }
151 double dist2_mm(Microcell* a, Microcell* b)
152 {
153  double x, y;
154  int l, m;
155 
156  l = (int)(a - Mcells);
157  m = (int)(b - Mcells);
158  if (P.DoUTM_coords)
159  return dist2UTM(P.in_microcells_.width_ * fabs((double)(l / P.get_number_of_micro_cells_high())), P.in_microcells_.height_ * fabs((double)(l % P.get_number_of_micro_cells_high())),
160  P.in_microcells_.width_ * fabs((double)(m / P.get_number_of_micro_cells_high())), P.in_microcells_.height_ * fabs((double)(m % P.get_number_of_micro_cells_high())));
161  else
162  {
163  x = P.in_microcells_.width_ * fabs((double)(l / P.get_number_of_micro_cells_high() - m / P.get_number_of_micro_cells_high()));
164  y = P.in_microcells_.height_ * fabs((double)(l % P.get_number_of_micro_cells_high() - m % P.get_number_of_micro_cells_high()));
165  if (P.DoPeriodicBoundaries)
166  {
167  if (x > P.in_degrees_.width_ * 0.5) x = P.in_degrees_.width_ - x;
168  if (y > P.in_degrees_.height_ * 0.5) y = P.in_degrees_.height_ - y;
169  }
170  return x * x + y * y;
171  }
172 }
173 
174 double dist2_raw(double ax, double ay, double bx, double by)
175 {
176  double x, y;
177 
178  if (P.DoUTM_coords)
179  return dist2UTM(ax, ay, bx, by);
180  else
181  {
182  x = fabs(ax - bx);
183  y = fabs(ay - by);
184  if (P.DoPeriodicBoundaries)
185  {
186  if (x > P.in_degrees_.width_ * 0.5) x = P.in_degrees_.width_ - x;
187  if (y > P.in_degrees_.height_ * 0.5) y = P.in_degrees_.height_ - y;
188  }
189  return x * x + y * y;
190  }
191 }
DomainSize in_degrees_
Size of spatial domain in degrees.
Definition: Param.h:106
The basic unit of the simulation and is associated to a geographical location.
Definition: Model.h:265
Definition: Model.h:15
DomainSize in_microcells_
Size of spatial domain in microcells.
Definition: Param.h:108
double height_
The height.
Definition: Param.h:24
Holds microcells.
Definition: Model.h:292
double width_
The width.
Definition: Param.h:21
DomainSize in_cells_
Size of spatial domain in cells.
Definition: Param.h:107