%% Visualize Antenna Field Strength Map on Earth
% This example shows how to calculate an antenna's field strength on flat
% earth and display it on the web map. The example also shows how to export
% the contour lines to KML, which can then be visualized on Earth browsers
% such as Google Earth(TM) and Google Maps(TM). This example requires the
% following product:
%
% * Mapping Toolbox(TM)

%   Copyright 2014-2016 The MathWorks, Inc.
%% Create Helix
% Use the default helix, and rotate it by 90 degrees so that most of its
% field is along the x-y plane i.e. on the earth surface.
f0   = 2.4e9;
ant  = helix('Tilt', -90);
show(ant);
view(-220, 30);

%% Calculate the Electric Field at Various Points in Space
% Assume the origin of the tilted antenna is located 15 m above the ground.
% Calculate the electric-field in a rectangular region covering 500 square
% km.
z = -15;
x = (-250:4:250)*1e3;
y = (-100:4:400)*1e3;

[X,Y]       = meshgrid(x,y);
TotalPoints = length(x)*length(y);
Points      = [X(:)'; Y(:)'; z*ones(1,numel(X))];

E = EHfields(ant, f0, Points); % Units: V/m

%% Calculate the Magnitude of Electric Field
% Compute the magnitude of the electric field.
Emag = zeros(1,TotalPoints);
for m=1:TotalPoints
    Emag(m) = norm(E(:,m));
end
Emag = 20*log10(reshape(Emag,length(y),length(x))); % Units: dBV/m
Emag = Emag + 120; % Units: dBuV/m

%% Plot the Electric Field
% Plot the magnitude of the electric field as a function of the x and y
% distance.
d_min = min(Emag(:));
d_max = max(Emag(:));
del   = (d_max-d_min)/12;
d_vec = round((d_min:del:d_max));

figure;
contourf(X*1e-3,Y*1e-3,Emag,d_vec,'showtext','on');
title('Field Strength (dB\muV/m) on flat Earth (1V tx)');
xlabel('lateral (km)'); 
ylabel('boresight (km)');
c = colorbar; 
set(get(c,'title'),'string','dB\muV/m');

%% Specify Antenna System Location and Orientation
% The first step in moving to an absolute coordinate system is to specify
% the location and orientation of the antenna. Assign a location of (42 N,
% 73 W) and a nominal height of 0 in WGS 84 geodetic coordinates.
lat0 =  42;
lon0 = -73;
h0 = 0;

%%
% Orient the antenna toward the southwest at an angle of -150 degrees
% (or, equivalently, +210 degrees) measured counterclockwise from East.
az = -150;

%%
% Because the antenna points along the local Y-axis (not the X-axis),
% subtract 90 degrees to determine the rotation of the local system with
% respect to xEast-yNorth. (This is the angle to the local X-axis from
% xEast measured counterclockwise.)
xyrot = wrapTo180(az - 90);

%% Contour Locally-Defined Input with Geodetic Output
% In this step, all the contouring happens at once. The results are a
% contourLines geoshape and a contourPolygons geoshape. Each feature
% (element) in contourLines consists of all the contour lines for a given
% level (identified by the feature property 'ContourLevel'). Likewise, each
% feature in contourPolygons is a multipart polygon enclosing the entire
% area filled for a given contour interval. The limits of the interval for
% a given polygon are stored in the feature properties LowerContourLevel
% and UpperContourLevel. Scale such that X and Y are (nominally) in meters.
% The default spheroid (which is the WGS 84 reference ellipsoid) will be in
% meters by default.

levelList = d_vec(d_vec >= 12); % Select contours that are complete
[contourLines, contourPolygons] = geocontourxy( ...
    X,Y,Emag,lat0,lon0,h0,'LevelList',levelList,'XYRotation',xyrot);
contourPolygons = contourPolygons(2:end); % Remove base shape
 
%% Display Geodetic Contour Lines in a Map Projection

% Set up base map.
[latlim, lonlim] = geoquadline( ...
    contourPolygons.Latitude, contourPolygons.Longitude);

latlim = latlim + [-1 1];
lonlim = lonlim + [-1 1];

states = shaperead('usastatelo.shp', ...
    'BoundingBox',[lonlim' latlim'],'UseGeoCoords',true);

figure
usamap(latlim, lonlim)
geoshow(states,'FaceColor',[0.8 0.8 0.8])

% Overlay contour lines.
minLevel = contourLines(1).ContourLevel;
maxLevel = contourLines(end).ContourLevel;
lineStyle = makesymbolspec('Line', ...
    {'ContourLevel',[minLevel maxLevel],'Color',parula});

h = geoshow(contourLines,'SymbolSpec',lineStyle);

% Add a marker to show the antenna location.
p = geopoint(lat0,lon0);
p.Description = 'System Location';
geoshow(p,'Marker','d','MarkerSize',10,'LineWidth',2, ...
    'MarkerEdgeColor','k','MarkerFaceColor','g')

% Add legend and title.
c = h.Children;
for k = 1:length(c)
    c(k).DisplayName = num2str(contourLines(k).ContourLevel);
end
l = legend(c, 'Location', 'BestOutside');
l.Title.String = 'Emag (dB\muV/m)';
title('Contours Projected on Map')

%% Display Geodetic Contour Polygons in Map Projection

% Set up base map.
figure
usamap(latlim, lonlim)
geoshow(states,'FaceColor',[0.8 0.8 0.8])

% Overlay contour polygons and antenna location.
cmap = parula(length(contourPolygons));
for k = 1:length(contourPolygons)
    lat = contourPolygons(k).Latitude;
    lon = contourPolygons(k).Longitude;
    geoshow(lat,lon,'Display','polygon', ...
        'FaceColor',cmap(k,:),'FaceAlpha',0.3,'EdgeColor','none')
end

% Overlay contour lines.
geoshow(contourLines.Latitude,contourLines.Longitude,'Color','black')

% Show system location.
geoshow(p,'Marker','d','MarkerSize',10,'LineWidth',2, ...
    'MarkerEdgeColor','k','MarkerFaceColor','g')

% Add title.
title('Filled Contours Projected on Map')

%% Display Contour Lines on Web Map
% The contour lines can be plotted on a webmap by using the following
% code. Bring up webmap and overlap the contour lines by using the function
% wmline and place a marker at the location by using wmmarker.
%%
%   webmap('world topographic map')
%   cmap = parula(length(contourLines));
%   wmline(contourLines,'Overlayname','Antenna Power 1', ...
%       'FeatureName','Antenna Contours',  ...
%       'Description','Antenna Contour','Color',cmap)
%   wmmarker(p,'Description','System Location','Overlayname','System Location')
%%
% 
% <<../webmap_lines.png>>
% 
%% Display Contour Polygons on Web Map
% The contour polygons can be plotted on a webmap by using the following
% code.
%%
%   cmap = parula(length(contourPolygons));
%   wmpolygon(contourPolygons,'Overlayname','Antenna Power 1', ...
%       'FeatureName','Antenna Contour Polygons',  ...
%       'Description','Antenna Contour Polygon', ...
%       'FaceColor',cmap,'FaceAlpha',0.3);
%   wmmarker(p,'Description','System Location','Overlayname','System Location')
%%
% 
% <<../webmap_polygon.png>>
% 
%% Export Contour Polygons to KML
% Using the command below, you can export the contour polygons to a KML
% file. This file can be used to visulaize the antenna pattern in Google
% Map(TM) or Google Earth(TM)
%%
%   kmlwrite('contourPolygons',contourPolygons(2:end),'AltitudeMode','clamp', ...
%   'Color',parula(length(contourPolygons)-1),'Alpha',0.3,'CutPolygon',false);
%
%%
