|
53 | 53 | % 27/09/2015 DB Separate the DEM loading, ask for number of domains. |
54 | 54 | % 28/09/2015 DB Also save the delay computation for a grid or z to include Pirate compartibility. |
55 | 55 | % 28/09/2015 DB Include multi-core from matlab |
| 56 | +% 07/07/2016 DB Redefine hydrostatic delay to be based on surface pressure. |
56 | 57 |
|
57 | 58 | save_complete=0; % save support information when 0 |
58 | 59 |
|
|
281 | 282 | ylist = ylist(ix); |
282 | 283 | latlist = latlist(ix); |
283 | 284 | lonlist = lonlist(ix); |
284 | | - |
285 | 285 | latlist=double(latlist); |
286 | 286 | lonlist=double(lonlist); |
287 | 287 |
|
| 288 | + numy = length(unique(latlist)); |
| 289 | + numx = length(unique(lonlist)); |
| 290 | + ulatlist = unique(latlist); |
| 291 | + ulonlist = unique(lonlist); |
| 292 | + uxlist = unique(xlist); |
| 293 | + uylist = unique(ylist); |
| 294 | + |
| 295 | + |
| 296 | + |
| 297 | + |
| 298 | + |
288 | 299 | if save_complete==0 |
289 | 300 | % saving the information for support plotting |
290 | 301 | wrf.wrf_lonlat = [lonlist latlist]; |
|
346 | 357 |
|
347 | 358 | %% Calculate Geometric Height, Z |
348 | 359 | Z = (H.*Re)./(g/g0.*Re - H); |
| 360 | + |
| 361 | + % Find middle of scene to work out glocal and Rlocal for use later |
| 362 | + midx = round(mean(uxlist)); |
| 363 | + midy = round(mean(uylist)); |
| 364 | + glocal = g(midy,midx,1); |
| 365 | + Rlocal = Re(midy,midx,1); |
| 366 | + |
| 367 | + |
349 | 368 |
|
350 | 369 | %% Find middle of scene to work out glocal and Rlocal for use later |
351 | 370 | cdslices = maxdem/vertres +1; |
|
354 | 373 | cdstack_wet = zeros(size(lonlist,1),cdslices); |
355 | 374 |
|
356 | 375 | XI=(0:zincr:zref)'; |
357 | | - %gh = glocal.*(Rlocal./(Rlocal+XI)).^2; %gravity with height for XI height range |
| 376 | + gh = glocal.*(Rlocal./(Rlocal+XI)).^2; %gravity with height for XI height range |
358 | 377 |
|
359 | 378 | % Interpolate Temp P and e from 0:20:15000 m |
360 | 379 | % then integrate using trapz to estimate delay as function of height |
|
383 | 402 | Lw = (10^-6).*-1*flipud(cumtrapz(flipud(XI),flipud(tmp1))); |
384 | 403 | % L is zenith delay one way in metres |
385 | 404 | tmp2 = k1.*YPI./YTI; |
386 | | - Ld = (10^-6).*-1*flipud(cumtrapz(flipud(XI),flipud(tmp2))); |
387 | | - |
| 405 | + %Ld = (10^-6).*-1*flipud(cumtrapz(flipud(XI),flipud(tmp2))); % This is using P/T expression (Hanssen, 2001) |
| 406 | + gm = glocal; |
| 407 | + Ld = (10^-6).*((k1*Rd/gm).*(YPI - YPI(zref/zincr +1))); % This is P0 expression (Hanssen, 2001) |
388 | 408 |
|
389 | 409 | % Interpolate important part (i.e. total delay at elevations |
390 | 410 | % less than maxdem) at high res i.e. vertres, and put in cdstack. |
|
0 commit comments