function [ ] = FootprintAreaRoll( h_km, lat_deg, lon_deg, az_deg, el_deg, azwidth_deg, elwidth_deg, roll_deg, kcolor )
                              
path( path, 'C:\Documents and Settings\guralp\My Documents\gural\MeteorWork\AnalysisMfiles' );

deg2rad = pi / 180;

%------ Site Info

rho = (pi/180) * 0.5 * sqrt( azwidth_deg^2 + elwidth_deg^2 );

phi = atan2( elwidth_deg, azwidth_deg );

[ ra, dec ] = AzimElev2RaDec( (az_deg) * deg2rad, el_deg * deg2rad, lat_deg * deg2rad, lon_deg * deg2rad );

[ razen, deczen ] = AzimElev2RaDec( 0.0 * deg2rad, 89.999 * deg2rad, lat_deg * deg2rad, lon_deg * deg2rad );

rhat(1) = cos(dec) * cos(ra);
rhat(2) = cos(dec) * sin(ra);
rhat(3) = sin(dec);

zhat(1) = cos(deczen) * cos(razen);
zhat(2) = cos(deczen) * sin(razen);
zhat(3) = sin(deczen);

azhat = cross( rhat, zhat );
azhat = azhat / norm(azhat);

elhat = cross( azhat, rhat );
elhat = elhat / norm(elhat);


%---- get first corner unit vector and resultant pointing vector

theta = phi + roll_deg * deg2rad;

aphat = azhat * cos(theta) + elhat * sin(theta);
aphat = aphat / norm(aphat);

pt = rhat * cos(rho) + aphat * sin(rho);
rhata = pt / norm(pt);

%---- get second corner unit vector and resultant pointing vector

theta = pi - phi + roll_deg * deg2rad;

aphat = azhat * cos(theta) + elhat * sin(theta);
aphat = aphat / norm(aphat);

pt = rhat * cos(rho) + aphat * sin(rho);
rhatb = pt / norm(pt);

%---- get first corner unit vector and resultant pointing vector

theta = pi + phi + roll_deg * deg2rad;

aphat = azhat * cos(theta) + elhat * sin(theta);
aphat = aphat / norm(aphat);

pt = rhat * cos(rho) + aphat * sin(rho);
rhatc = pt / norm(pt);

%---- get first corner unit vector and resultant pointing vector

theta = -phi + roll_deg * deg2rad;

aphat = azhat * cos(theta) + elhat * sin(theta);
aphat = aphat / norm(aphat);

pt = rhat * cos(rho) + aphat * sin(rho);
rhatd = pt / norm(pt);



[ x_km, y_km, z_km ] = LLA2ECEF( lat_deg, lon_deg, h_km );

site1(1) = x_km;
site1(2) = y_km;
site1(3) = z_km;


%------ plot the areas that intesect at h kilometers

h = 90;


xlabel('Longitude (deg)');
ylabel('Latitude (deg)');

for range=1:0.1:1000
    r = site1 + range * rhata;
    [ lat1_deg, lon1_deg, alt_km ] = ECEF2LLA( r(1), r(2), r(3) );
    if alt_km >= h
        break;
    end
end
lat(1) = lat1_deg;
lon(1) = lon1_deg;
    
for range=1:0.1:1000
    r = site1 + range * rhatb;
    [ lat1_deg, lon1_deg, alt_km ] = ECEF2LLA( r(1), r(2), r(3) );
    if alt_km >= h
        break;
    end
end
lat(2) = lat1_deg;
lon(2) = lon1_deg;
    
for range=1:0.1:1000
    r = site1 + range * rhatc;
    [ lat1_deg, lon1_deg, alt_km ] = ECEF2LLA( r(1), r(2), r(3) );
    if alt_km >= h
        break;
    end
end
lat(3) = lat1_deg;
lon(3) = lon1_deg;
    
for range=1:0.1:1000
    r = site1 + range * rhatd;
    [ lat1_deg, lon1_deg, alt_km ] = ECEF2LLA( r(1), r(2), r(3) );
    if alt_km >= h
        break;
    end
end
lat(4) = lat1_deg;
lon(4) = lon1_deg;
    
for range=1:0.1:1000
    r = site1 + range * rhata;
    [ lat1_deg, lon1_deg, alt_km ] = ECEF2LLA( r(1), r(2), r(3) );
    if alt_km >= h
        break;
    end
end
lat(5) = lat1_deg;
lon(5) = lon1_deg;


if kcolor==1
    plot( lon, lat, 'k' );  
    plot( lon_deg, lat_deg, 'ko' );
elseif  kcolor==2
    plot( lon, lat, 'm' );  
    plot( lon_deg, lat_deg, 'mo' );
elseif  kcolor==3
    plot( lon, lat, 'b' );  
    plot( lon_deg, lat_deg, 'bo' );
elseif  kcolor==4
    plot( lon, lat, 'g' );  
    plot( lon_deg, lat_deg, 'go' );
elseif  kcolor==5
    plot( lon, lat, 'c' );  
    plot( lon_deg, lat_deg, 'co' );
elseif  kcolor==6
    plot( lon, lat, 'r' );  
    plot( lon_deg, lat_deg, 'ro' );
else
    plot( lon, lat, 'k' );  
    plot( lon_deg, lat_deg, 'ko' );
end


%=========== Output kml file format

% 		<Polygon>
% 			<outerBoundaryIs>
% 				<LinearRing>
% 				  <coordinates>
% 				   -78.8765,39.9835,0
% 				   -78.0476,39.4603,0
% 				   -77.6830,39.6617,0
% 				   -78.0626,40.4351,0
% 				   -78.8765,39.9835,0
% 				  </coordinates>
% 				</LinearRing>
% 			</outerBoundaryIs>
% 		</Polygon>


disp('      <Polygon>');
disp('      	<outerBoundaryIs>');
disp('      		<LinearRing>');
disp('      		  <coordinates>');
for kk=1:5
    s = sprintf('\t\t		   %f,%f,0',lon(kk),lat(kk));
    disp(s);
end
disp('      		  </coordinates>');
disp('      		</LinearRing>');
disp('      	</outerBoundaryIs>');
disp('      </Polygon>');
disp('      ');