clear all; close all; clc;
%%  Example 1: Gaussian mountain (theoretical)

hm = 10000;          % mnt height
xm = 300000;         % mountain location in domain 
ym = 300000;
am = 150000;         % mountain half-width (m)
dx = 500;            % grid spacing

[x, y] = meshgrid(0:dx:(xm*2), 0:dx:(xm*2));
zs = hm*exp(-(((x-xm)./am).^2 + ((y-ym)./am).^2));  % Gaussian
P = LT3d(x, y, zs, 15, 0, 'hw', 1000);

% Plot
figure()
imagesc(P)
hold on
contour(zs, 10, 'w')
caxis([0 30])
colorbar;
set(gca, 'ydir', 'normal');

%% Example 2: Realworld example with terrain data for the Nauelbuta mountains

zz = srtmread([-74 -72], [-38.5 -36.5]);
[xx,yy] = meshgrid((0:(size(zz,2)-1))*100, (0:(size(zz,2)-1))*100);
lon = linspace(-74,-72, size(zz,2));
lat = linspace(-38.5,-36.5, size(zz,2));
[LON, LAT] = meshgrid(lon, lat);

P = LT3d(xx, yy, zz, 10, -5, 'hw', 200);

% Plot
figure()
imagesc(lon,lat,P)
hold on
contour(LON,LAT,zz,[1,500,1000,1500],'w')
caxis([0 10])
colorbar;
set(gca, 'ydir', 'normal');

