← Back to branch summary

~fluidity-core/fluidity/sea-ice-coupling

« back to all changes in this revision

Viewing changes to bathymetry/import_bath_data.cpp

  • Committer: Simon Mouradian
  • Date: 2012-05-21 21:54:54 UTC
  • mfrom: (3520.32.268 fluidity)
  • Revision ID: mouradian@gmail.com-20120521215454-xhzgwfcb0bnzuww1
merge from lp:fluidity

Show diffs side-by-side

added added

removed removed

Lines of Context:
bathymetry/import_bath_data.cpp
38
38
 
39
39
using namespace std;
40
40
 
41
 
void get_lldepth(double x, double y, double z, double &longitude, double &latitude, double &depth){
 
41
void get_ll(double x, double y, double z, double &longitude, double &latitude){
42
42
 
43
43
  const double pi=4.0*atanl(1.0);
44
 
  const double earth_radius=6378100.0;
45
44
  const double rad_to_deg=180.0/pi;
46
45
  
47
46
  double r = sqrt(x*x+y*y+z*z);
48
47
  longitude = atan2(y, x);
49
48
  latitude = acos(z/r);
50
 
  depth = earth_radius - r;
51
49
  longitude*=rad_to_deg;
52
50
  latitude = 90.0 - latitude*rad_to_deg;
53
51
}
75
73
          x[i] = X[i];
76
74
          y[i] = Y[i];
77
75
          z[i] = Z[i];
78
 
          double longitude, latitude, depth_tmp;
79
 
          get_lldepth(x[i], y[i], z[i], longitude, latitude, depth_tmp);
 
76
          double longitude, latitude;
 
77
          get_ll(x[i], y[i], z[i], longitude, latitude);
80
78
          if(map.HasPoint(longitude, latitude)){
81
79
                height[i]=map.GetValue(longitude, latitude);
82
80
          }else{
127
125
        }
128
126
}
129
127
 
130
 
#ifdef DEPTH_UNITTEST
131
 
// This is a simple test program to check the depth returned for three given lon,lat locations
132
 
int main(int argc, char **argv){
133
 
 
134
 
  double pi=4.0*atanl(1.0);
135
 
  double deg_to_rad=pi/180.0;
136
 
  double depth[3];
137
 
  int columns=3;
138
 
 
139
 
  string data_file="/data/maps/gridone.grd";
140
 
 
141
 
  double xloc[3], yloc[3], zloc[3], longitude[3], latitude[3];
142
 
 
143
 
  double r=6378100.0; // Approx. radius of Earth
144
 
 
145
 
  // 1st location
146
 
  longitude[0]=5.974731;
147
 
  latitude[0]=38.387311;
148
 
 
149
 
  xloc[0]=r*cos(latitude[0]*deg_to_rad)*cos(longitude[0]*deg_to_rad);
150
 
  yloc[0]=r*cos(latitude[0]*deg_to_rad)*sin(longitude[0]*deg_to_rad);
151
 
  zloc[0]=r*sin(latitude[0]*deg_to_rad);
152
 
 
153
 
  // 2nd location
154
 
  longitude[1]=6.974731;
155
 
  latitude[1]=39.387311;
156
 
 
157
 
  xloc[1]=r*cos(latitude[1]*deg_to_rad)*cos(longitude[1]*deg_to_rad);
158
 
  yloc[1]=r*cos(latitude[1]*deg_to_rad)*sin(longitude[1]*deg_to_rad);
159
 
  zloc[1]=r*sin(latitude[1]*deg_to_rad);
160
 
 
161
 
  // 3rd location
162
 
  longitude[2]=5.974731;
163
 
  latitude[2]=40.387311;
164
 
 
165
 
  xloc[2]=r*cos(latitude[2]*deg_to_rad)*cos(longitude[2]*deg_to_rad);
166
 
  yloc[2]=r*cos(latitude[2]*deg_to_rad)*sin(longitude[2]*deg_to_rad);
167
 
  zloc[2]=r*sin(latitude[2]*deg_to_rad);
168
 
 
169
 
  cout << "Retrieving bathymetry data.\n";
170
 
  set_from_map_fc(data_file.c_str(), xloc, yloc, zloc, depth, &columns);
171
 
 
172
 
  for (int i = 0; i < columns; i++) {
173
 
    cout << "The depth at ( " << longitude[i] << " , " << latitude[i] << " ) is " << depth[i] << endl;
174
 
  }
175
 
  cout << "End of unit test.\n";  
176
 
  
177
 
}
178
 
#endif
179
 
 
180
 
#ifdef DEPTH_UNITTEST2
181
 
// This is a simple test program to check the depth returned for three given x,y beta plane locations
182
 
int main(int argc, char **argv){
183
 
 
184
 
  double depth[3];
185
 
  int columns=3;
186
 
 
187
 
  string data_file="/data/myfluidity/tests/MonaiValley/outfileMonaiValley.grd";
188
 
 
189
 
  double xloc[3], yloc[3];
190
 
  
191
 
  // The Monaia Valley data bounds are
192
 
  // x:actual_range = 0., 5.4916666667
193
 
  // y:actual_range = 0., 3.4083333333
194
 
 
195
 
  // 1st location
196
 
  xloc[0]=1.0;
197
 
  yloc[0]=1.0;
198
 
 
199
 
  // 2nd location
200
 
  xloc[1]=2.0;
201
 
  yloc[1]=2.0;
202
 
 
203
 
  // 3rd location
204
 
  xloc[2]=3.0;
205
 
  yloc[2]=3.0;
206
 
 
207
 
  cout << "Retrieving bathymetry data.\n";
208
 
  set_from_map_beta_fc(data_file.c_str(), xloc, yloc, depth, &columns);
209
 
 
210
 
  for (int i = 0; i < columns; i++) {
211
 
    cout << "The depth at ( " << xloc[i] << " , " << yloc[i] << " ) is " << depth[i] << endl;
212
 
  }
213
 
  cout << "End of unit test.\n";  
214
 
  
215
 
}
216
 
#endif
217