% Directed Energy Power Density Levels
% ------------------------------------

clear;clc;clf;

% Enter System Parameters

D=3;				   % Antenna Diameter - m
Pd=1e5;			   % Near Field Power Density - w/cm^2
R1=10;   		   % Range 1 - km
r=[1 4];			   % Rain Rate - mm/hr

% Enter Frequency Range

f=1:30;		   	% Frequency - GHz

% Compute Power in Near Field

Pnf=Pd*pi*(3*100)^2/4;

% Compute Antenna Gain

wl=3e8./(f*1e9);  % Wavelength - m
G=(pi*D./wl).^2;

% Compute Power Density in Far Field

P=Pnf*G/(4*pi*(1000*R1)^2);  

% Include Effects of Atmosphere 

% Frequency Points - GHz

fx=[.4 1.3 3 5.5 10 15 22 35];

% Attentuation ( Two-Way ) - dB/km

att=[.01 .012 .015 .017 .024 .055 .30 .14];

fac=R1*interp1(fx,att,f,'spline');

Px=P.*10.^(-5*fac);

% Plot Power Density

semilogy(f,P,f,Px);grid;
xlabel('Frequency - GHz');
ylabel('Power Density -  w /  m  ^2');
title(['Maximum Power Density  ( R =',num2str(R1),' km )']);
axis([ 0 30 1e-60 1e10]);
hold on;

% Include Effects of Percipitation

att1=[0 .0003 .0013 .008 .037 .083 .23 .57];

fac1=R1*interp1(fx,att1,f,'spline');

Px1=Px.*10.^(-5*fac1*r(1));
Px2=Px.*10.^(-5*fac1*r(2));
semilogy(f,Px1,f,Px2);