%script to generate an ARC shape file from SBDR info structure
function sbdr_make_shapefile_v2(filename);

%define small major axis threshold (degrees);
small_thresh = 10;
scnt = 1;

info=-1;
if length(strfind(filename,'.mat')) > 0;    
    %load the .mat file
    load(filename);
else;
    disp('UnRecognized Data File Format!');
    disp(['Filename: ' filename]);
    return;
end;

%check the format of the info structure
field_names = fields(info);
if ~strcmp(field_names{1},lower(field_names{1}));
    for i=1:length(field_names);
        eval(['info2.' lower( field_names{i} ) ...
            ' = info.' field_names{i} ';']);
        eval(['info = info2;']);
        clear info2;
        save(filename,'info');
    end;
end;
        

%create output file names
file_out = [filename( 1:(strfind( filename, '.mat')-1)  ) '_Rad.shp'];
file_out2 = [filename( 1:(strfind( filename, '.mat')-1)  ) '_Rad_centroid.shp'];
file_out3 = [filename( 1:(strfind( filename, '.mat')-1)  ) '_Rad_center.shp'];
file_out4 = [filename( 1:(strfind( filename, '.mat')-1)  ) '_Rad_small.shp'];
file2_out = [filename( 1:(strfind( filename, '.mat')-1)  ) '_Rad.shp'];
file2_out2 = [filename( 1:(strfind( filename, '.mat')-1)  ) '_Scatt_centroid.shp'];
file2_out3 = [filename( 1:(strfind( filename, '.mat')-1)  ) '_Scatt_center.shp'];
file2_out4 = [filename( 1:(strfind( filename, '.mat')-1)  ) '_Scatt_small.shp'];

%figure out the number of good points in the recond
science_qual_flag = (dec2bin( info.science_qual_flag ));
if size(science_qual_flag,2) ~= 10;
    disp('Bad Science Qual Flag: Using All of the Available Data!');
    ind = find( ( info.antenna_temp ~= 0 ) & ...
                ( (info.pass_ellipse_pt1_lat + ...
                   info.pass_ellipse_pt1_lon + ...
                   info.pass_ellipse_pt2_lat + ...
                   info.pass_ellipse_pt2_lon + ...
                   info.pass_ellipse_pt3_lat + ...
                   info.pass_ellipse_pt3_lon + ...
                   info.pass_ellipse_pt4_lat + ...
                   info.pass_ellipse_pt4_lon) ~= 0 ) );
    ind2 = find( ( info.sigma0_corrected ~= 0 ) & ...
                ( (info.act_ellipse_pt1_lat + ...
                   info.act_ellipse_pt1_lon + ...
                   info.act_ellipse_pt2_lat + ...
                   info.act_ellipse_pt2_lon + ...
                   info.act_ellipse_pt3_lat + ...
                   info.act_ellipse_pt3_lon + ...
                   info.act_ellipse_pt4_lat + ...
                   info.act_ellipse_pt4_lon) ~= 0 ) );
else;    
    %Radiometry
    ind = find( (science_qual_flag(:,10) == '0') & ...
                (science_qual_flag(:,6) == '0') & ...
                (science_qual_flag(:,5) == '0') & ...
                (science_qual_flag(:,4) == '0') & ...
                ( (info.pass_ellipse_pt1_lat + ...
                   info.pass_ellipse_pt1_lon + ...
                   info.pass_ellipse_pt2_lat + ...
                   info.pass_ellipse_pt2_lon + ...
                   info.pass_ellipse_pt3_lat + ...
                   info.pass_ellipse_pt3_lon + ...
                   info.pass_ellipse_pt4_lat + ...
                   info.pass_ellipse_pt4_lon) ~= 0 ) );
       %Scatterometry
       ind2 = find( (science_qual_flag(:,9) == '0') & ...
           (science_qual_flag(:,7) == '0') & ...
           (science_qual_flag(:,3) == '0') & ...
           (science_qual_flag(:,2) == '0') & ...
           ( (info.act_ellipse_pt1_lat + ...
           info.act_ellipse_pt1_lon + ...
           info.act_ellipse_pt2_lat + ...
           info.act_ellipse_pt2_lon + ...
           info.act_ellipse_pt3_lat + ...
           info.act_ellipse_pt3_lon + ...
           info.act_ellipse_pt4_lat + ...
           info.act_ellipse_pt4_lon) ~= 0 ) );             
end;

if length(ind) == 0;
    disp(['Filename: ' filename]);
    disp('Error: No Useful Radiometer Data!');
	badrad = 1;
else;
    num_record = length(ind);
    disp(['Filename: ' filename]);
    disp(['Creating Shape File With ' num2str( num_record ) ' Rad. Entries']);
    badrad = 0;
end;

if length(ind2) == 0;
    disp('Error: No Useful Scatterometry Data!');
    badscatt = 1;   
else;
    num_record2 = length(ind2);
    disp(['Creating Shape File With ' num2str( num_record ) ' Scat. Entries']);
    badscatt = 1;
end;
if ( badrad & badscatt );
    return;
end;

if badrad == 0;
%step through each footprint and creat 100 lat lon points around the
%ellipse
for i=1:num_record;
    %get the ellipse center and semi-major axis, eccentricity, and
    %azimuth
    pt1x = info.pass_ellipse_pt1_lat(ind(i));
    pt1y = info.pass_ellipse_pt1_lon(ind(i));
    pt2x = info.pass_ellipse_pt2_lat(ind(i));
    pt2y = info.pass_ellipse_pt2_lon(ind(i));
    pt3x = info.pass_ellipse_pt3_lat(ind(i));
    pt3y = info.pass_ellipse_pt3_lon(ind(i));
    pt4x = info.pass_ellipse_pt4_lat(ind(i));
    pt4y = info.pass_ellipse_pt4_lon(ind(i)); 
    
    %convert lons into -180 to 180
    if pt1y > 180; pt1y = pt1y - 360; end;
    if pt2y > 180; pt2y = pt2y - 360; end;
    if pt3y > 180; pt3y = pt3y - 360; end;
    if pt4y > 180; pt4y = pt4y - 360; end;
    if pt1y < -180; pt1y = pt1y + 360; end;
    if pt2y < -180; pt2y = pt2y + 360; end;
    if pt3y < -180; pt3y = pt3y + 360; end;
    if pt4y < -180; pt4y = pt4y + 360; end;
    
    ptx = [pt1x,pt2x,pt3x,pt4x];
    pty = -1.*[pt1y,pt2y,pt3y,pt4y];
    
    %%see if we are cross 180 or 360
    if ( min( abs( pty - 360 ) ) < 2 ) | (min(abs(pty)) < 2);
        %convert to -180 to 180
        pty(pty > 180) = pty(pty > 180) -360;
    elseif ( min( abs( pty -180) ) < 2 );
        %convert to 0 to 360
        pty(pty < 0) = pty(pty < 0) + 360;
    end;

    %make sure these are single records
    pt1y = pty(1);
    pt2y = pty(2);
    pt3y = pty(3);
    pt4y = pty(4);
    pt1x = ptx(1);
    pt2x = ptx(2);
    pt3x = ptx(3);
    pt4x = ptx(4);
        
    %find intersection point using triangulation
    xvals = cos( [pt1y,pt3y,pt2y,pt4y] .*pi/180) .* ...
        cos( [pt1x,pt3x,pt2x,pt4x] .*pi/180);
    yvals = sin( [pt1y,pt3y,pt2y,pt4y] .*pi/180) .* ...
        cos( [pt1x,pt3x,pt2x,pt4x] .*pi/180);
    zvals = sin( [pt1x,pt3x,pt2x,pt4x] .*pi/180);
    A = [xvals(1),yvals(1),zvals(1)];
    B  = [xvals(2),yvals(2),zvals(2)];
    C  = [xvals(3),yvals(3),zvals(3)];
    D =  [xvals(4),yvals(4),zvals(4)];
    N = cross(B-A,C-A);
    N = N./norm(N);
    latc1 = -asin(N(3)).*180/pi;
    lonc1= atan2(N(2),N(1))*180/pi-180;
    N2 = cross(B-A,D-A);
    N2 = N2./norm(N2);
    latc2 = -asin(N2(3))*180/pi;
    lonc2 = atan2(N2(2),N2(1))*180/pi-180;
    N3 = cross(C-B,D-B);
    N3 = N3./norm(N3);
    latc3 = -asin(N3(3)).*180/pi;
    lonc3= atan2(N3(2),N3(1))*180/pi-180;
    N4 = cross(D-C,A-C);
    N4 = N4./norm(N4);
    latc4 = -asin(N4(3))*180/pi;
    lonc4 = atan2(N4(2),N4(1))*180/pi-180;
    lat0 = mean([latc1,latc2,latc3,latc4]);
    lon0 = mean([lonc1,lonc2,lonc3,lonc4]);
    if lon0 < (-180); lon0 = lon0 + 360; end;
    if lon0 > ( 180); lon0 = lon0 - 360; end;
    lat0s(i) = lat0;
    lon0s(i) = lon0;
    latcens(i) = info.pass_centroid_lat(ind(i));
    loncens(i) = -info.pass_centroid_lon(ind(i));
    if loncens(i) < (-180); loncens(i) = loncens(i) + 360; end;
    if loncens(i) > ( 180); loncens(i) = loncens(i) - 360; end;
    
    %derive semi-major and semi-minor axis for geographic ellipse
    a = sqrt( (pt1x-pt2x).^2 + (pt1y-pt2y).^2 );
    b = sqrt( (pt3x-pt4x).^2 + (pt3y-pt4y).^2 );
    c = sqrt( (pt1x-pt3x).^2 + (pt1y-pt3y).^2 );
    d = sqrt( (pt1x-pt4x).^2 + (pt1y-pt4y).^2 );
    e = sqrt( (pt2x-pt3x).^2 + (pt2y-pt3y).^2 );
    f = sqrt( (pt2x-pt4x).^2 + (pt2y-pt4y).^2 );
    [major,mind] = max([a,b,c,d,e,f]);
    switch mind
        case 1
            minor = b;
            az2 = atan( (pt1y-pt2y)./(pt1x-pt2x) );
        case 2
            minor = a;
            az2 = atan( (pt4y-pt3y)./(pt4x-pt3x) );
        case 3
            minor = f;
            az2 = atan( (pt1y-pt3y)./(pt1x-pt3x) );
        case 4
            minor = e;
            az2 = atan( (pt4y-pt1y)./(pt4x-pt1x) );
        case 5
            minor = d;
            az2 = atan( (pt2y-pt3y)./(pt2x-pt3x) );
        case 6;
            minor = c;
            az2 = atan( (pt4y-pt2y)./(pt4x-pt2x) );
    end;       
    ecc = sqrt( 1 - (minor./major).^2);

    if major < small_thresh;
        yes_small = 1;
    else;
        yes_small = 0;
    end;
    
    %%find angular offsets
    %         %normalize semimajor axis
    %         n1 =[ pt1x - pt2x, pt1y-pt2y ]./sqrt( (pt1x-pt2x)^2+(pt1y-pt2y)^2);
    %         off1 = acos( n1(2) ) * 180/pi;
    %         %normalize semiminor axis
    %         n2 =[ pt3x - pt4x, pt3y-pt4y ]./sqrt( (pt3x-pt4x)^2+(pt3y-pt4y)^2);
    %         off2 = acos( n2(2) ) * 180/pi;
    
    %get the reported azimuthal look direction
    az = (-1).*(info.pass_azimuth_angle(ind(i)));
    
    %create the ellipse    
    t = linspace(-pi,pi,100);
    lat = lat0 + major./2.*cos(t).*cos(az2) - ...
                 minor./2.*sin(t).*sin(az2);             
    lon = lon0 + major./2.*cos(t).*sin(az2) + ...
                 minor./2.*sin(t).*cos(az2);
                 
%     %create the ellipse
%     [lat,lon] = ellipse1(lat0,lon0,[major./2,ecc],az);

    %convert the longitude values into -180 to 180 
    lon(lon > 180) = lon(lon > 180) - 360;
    lon0s(lon0s > 180) = lon0s(lon0s > 180) - 360;
    loncens(loncens > 180) = loncens(loncens > 180) - 360;
    
    %if we pass through 180 in the middle fix it
    if ( max(lon) > 170 ) & ( min(lon) < -170 );
        lon(lon < 0) = lon(lon < 0) + 360;
    end;
    

%     %fix the sign conventions
%     loncens = -loncens;
%     lat = -lat;
%     lon = (-lon)+180;
%     lon0s(i) = -(lon0s(i))+180;
%     lat0s(i) = -lat0s(i);
       
    %populate the structure file
    struct(i).BoundingBox = ...
        [min(lat) max(lat); min(lon) max(lon)];
    struct(i).Y = lat;
    struct(i).X = lon;
    struct(i).Geometry = 'Polygon';
    struct(i).Id = i;
    struct(i).Energy = info.total_echo_energy(ind(i))*1e16;
    struct(i).Energy_Units = 'Joules * 1e16';
    struct(i).Sigma0_Cor = info.sigma0_corrected(ind(i));
    struct(i).Sigma0_UnCor = info.sigma0_uncorrected(ind(i));
    struct(i).Num_Pulse = info.num_pulses_received(ind(i));
    struct(i).Gain = info.system_gain(ind(i));
    struct(i).Antenna = info.antenna_temp(ind(i));
    struct(i).Antenna_Std = info.ant_temp_std(ind(i));
    struct(i).Receiver = info.receiver_temp(ind(i));
    struct(i).Pass_Pol_Angle = info.pass_pol_angle(ind(i));
    struct(i).Pass_Azimuth_Angle = info.pass_azimuth_angle(ind(i));
    struct(i).Pass_Emission_Angle = info.pass_emission_angle(ind(i));
    struct(i).Act_Pol_Angle = info.act_pol_angle(ind(i));
    struct(i).Act_Incidence_Angle = info.act_indcidence_angle(ind(i));
    struct(i).Act_Azimuth_Angle = info.act_azimuth_angle(ind(i));
    struct(i).X_Factor = info.x_factor(ind(i));
    struct(i).Radar_Mode = info.radar_mode(ind(i));
    struct2(i).X = loncens(i);
    struct2(i).Y = latcens(i);
    struct2(i).Geometry = 'Point';
    struct2(i).Id = i;
    struct2(i).Energy = info.total_echo_energy(ind(i))*1e16;
    struct2(i).Energy_Units = 'Joules * 1e16';
    struct2(i).Sigma0_Cor = info.sigma0_corrected(ind(i));
    struct2(i).Sigma0_UnCor = info.sigma0_uncorrected(ind(i));
    struct2(i).Num_Pulse = info.num_pulses_received(ind(i));
    struct2(i).Gain = info.system_gain(ind(i));   
    struct2(i).Antenna = info.antenna_temp(ind(i));
    struct2(i).Antenna_Std = info.ant_temp_std(ind(i));
    struct2(i).Receiver = info.receiver_temp(ind(i));
    struct2(i).Pass_Pol_Angle = info.pass_pol_angle(ind(i));
    struct2(i).Pass_Azimuth_Angle = info.pass_azimuth_angle(ind(i));
    struct2(i).Pass_Emission_Angle = info.pass_emission_angle(ind(i));
    struct2(i).Act_Pol_Angle = info.act_pol_angle(ind(i));
    struct2(i).Act_Incidence_Angle = info.act_indcidence_angle(ind(i));
    struct2(i).Act_Azimuth_Angle = info.act_azimuth_angle(ind(i));
    struct2(i).X_Factor = info.x_factor(ind(i));
    struct2(i).Radar_Mode = info.radar_mode(ind(i));
    struct3(i).X = lon0s(i);
    struct3(i).Y = lat0s(i);
    struct3(i).Geometry = 'Point';
    struct3(i).Id = i;
    struct3(i).Energy = info.total_echo_energy(ind(i))*1e16;
    struct3(i).Energy_Units = 'Joules * 1e16';
    struct3(i).Sigma0_Cor = info.sigma0_corrected(ind(i));
    struct3(i).Sigma0_UnCor = info.sigma0_uncorrected(ind(i));
    struct3(i).Num_Pulse = info.num_pulses_received(ind(i));
    struct3(i).Gain = info.system_gain(ind(i));     
    struct3(i).Antenna = info.antenna_temp(ind(i));
    struct3(i).Antenna_Std = info.ant_temp_std(ind(i));
    struct3(i).Receiver = info.receiver_temp(ind(i));
    struct3(i).Pass_Pol_Angle = info.pass_pol_angle(ind(i));
    struct3(i).Pass_Azimuth_Angle = info.pass_azimuth_angle(ind(i));
    struct3(i).X_Factor = info.x_factor(ind(i));
    struct3(i).Pass_Emission_Angle = info.pass_emission_angle(ind(i));
    struct3(i).Act_Pol_Angle = info.act_pol_angle(ind(i));
    struct3(i).Act_Incidence_Angle = info.act_indcidence_angle(ind(i));
    struct3(i).Act_Azimuth_Angle = info.act_azimuth_angle(ind(i));
    struct3(i).Radar_Mode = info.radar_mode(ind(i));
    if yes_small;
        struct4(scnt).BoundingBox = ...
            [min(lat) max(lat); min(lon) max(lon)];
        struct4(scnt).Y = lat;
        struct4(scnt).X = lon;
        struct4(scnt).Geometry = 'Polygon';
        struct4(scnt).Id = i;
        struct4(scnt).Gain = info.system_gain(ind(i));        
        struct4(scnt).Energy = info.total_echo_energy(ind(i))*1e16;
        struct4(scnt).Energy_Units = 'Joules * 1e16';
        struct4(scnt).Sigma0_Cor = info.sigma0_corrected(ind(i));
        struct4(scnt).Sigma0_UnCor = info.sigma0_uncorrected(ind(i));
        struct4(scnt).Num_Pulse = info.num_pulses_received(ind(i));
        struct4(scnt).Gain = info.system_gain(ind(i));
        struct4(scnt).Antenna = info.antenna_temp(ind(i));
        struct4(scnt).Antenna_Std = info.ant_temp_std(ind(i));
        struct4(scnt).Receiver = info.receiver_temp(ind(i));
        struct4(scnt).Pass_Pol_Angle = info.pass_pol_angle(ind(i));
        struct4(scnt).Pass_Azimuth_Angle = info.pass_azimuth_angle(ind(i));
        struct4(scnt).X_Factor = info.x_factor(ind(i));
        struct4(scnt).Pass_Emission_Angle = info.pass_emission_angle(ind(i));
        struct4(scnt).Act_Pol_Angle = info.act_pol_angle(ind(i));
        struct4(scnt).Act_Incidence_Angle = info.act_indcidence_angle(ind(i));
        struct4(scnt).Act_Azimuth_Angle = info.act_azimuth_angle(ind(i));
        struct4(scnt).Radar_Mode = info.radar_mode(ind(i));
        scnt = scnt+1;
    end;
end;

%write out the shape file
shapewrite(struct,file_out);
shapewrite(struct2,file_out2);
shapewrite(struct3,file_out3);
shapewrite(struct4,file_out4);

%create projection files for the new shape files
%open and write a projection file for titan to go along with it
out_file = file_out;
pos = strfind( out_file, '.shp' );
prj_file = [out_file(1:pos) 'prj'];
fid = fopen(prj_file,'w');
string = ['GEOGCS["GCS_Titan_2000",DATUM["D_Titan_2000",'...
    'SPHEROID["Titan_2000_IAU_IAG",2575000.0,0.0]],PRIMEM["Reference_Meridian",0.0],',...
    'UNIT["Degree",0.0174532925199433]]' '\n'];
fprintf(fid,string);
fclose(fid);
out_file = file_out2;
pos = strfind( out_file, '.shp' );
prj_file = [out_file(1:pos) 'prj'];
fid = fopen(prj_file,'w');
string = ['GEOGCS["GCS_Titan_2000",DATUM["D_Titan_2000",'...
    'SPHEROID["Titan_2000_IAU_IAG",2575000.0,0.0]],PRIMEM["Reference_Meridian",0.0],',...
    'UNIT["Degree",0.0174532925199433]]' '\n'];
fprintf(fid,string);
fclose(fid);
out_file = file_out3;
pos = strfind( out_file, '.shp' );
prj_file = [out_file(1:pos) 'prj'];
fid = fopen(prj_file,'w');
string = ['GEOGCS["GCS_Titan_2000",DATUM["D_Titan_2000",'...
    'SPHEROID["Titan_2000_IAU_IAG",2575000.0,0.0]],PRIMEM["Reference_Meridian",0.0],',...
    'UNIT["Degree",0.0174532925199433]]' '\n'];
fprintf(fid,string);
fclose(fid);
out_file = file_out4;
pos = strfind( out_file, '.shp' );
prj_file = [out_file(1:pos) 'prj'];
fid = fopen(prj_file,'w');
string = ['GEOGCS["GCS_Titan_2000",DATUM["D_Titan_2000",'...
    'SPHEROID["Titan_2000_IAU_IAG",2575000.0,0.0]],PRIMEM["Reference_Meridian",0.0],',...
    'UNIT["Degree",0.0174532925199433]]' '\n'];
fprintf(fid,string);
fclose(fid);

disp(['Done with Radiometry For ' filename]);
end;

if badscatt == 0;
num_record = num_record2;
scnt = 1;
%step through each footprint and creat 100 lat lon points around the
%ellipse
for i=1:num_record;
    %get the ellipse center and semi-major axis, eccentricity, and
    %azimuth
    pt1x = info.act_ellipse_pt1_lat(ind(i));
    pt1y = info.act_ellipse_pt1_lon(ind(i));
    pt2x = info.act_ellipse_pt2_lat(ind(i));
    pt2y = info.act_ellipse_pt2_lon(ind(i));
    pt3x = info.act_ellipse_pt3_lat(ind(i));
    pt3y = info.act_ellipse_pt3_lon(ind(i));
    pt4x = info.act_ellipse_pt4_lat(ind(i));
    pt4y = info.act_ellipse_pt4_lon(ind(i)); 
    
    %convert lons into -180 to 180
    if pt1y > 180; pt1y = pt1y - 360; end;
    if pt2y > 180; pt2y = pt2y - 360; end;
    if pt3y > 180; pt3y = pt3y - 360; end;
    if pt4y > 180; pt4y = pt4y - 360; end;
    if pt1y < -180; pt1y = pt1y + 360; end;
    if pt2y < -180; pt2y = pt2y + 360; end;
    if pt3y < -180; pt3y = pt3y + 360; end;
    if pt4y < -180; pt4y = pt4y + 360; end;
    
    ptx = [pt1x,pt2x,pt3x,pt4x];
    pty = -1.*[pt1y,pt2y,pt3y,pt4y];
    
    %%see if we are cross 180 or 360
    if ( min( abs( pty - 360 ) ) < 2 ) | (min(abs(pty)) < 2);
        %convert to -180 to 180
        pty(pty > 180) = pty(pty > 180) -360;
    elseif ( min( abs( pty -180) ) < 2 );
        %convert to 0 to 360
        pty(pty < 0) = pty(pty < 0) + 360;
    end;

    %make sure these are single records
    pt1y = pty(1);
    pt2y = pty(2);
    pt3y = pty(3);
    pt4y = pty(4);
    pt1x = ptx(1);
    pt2x = ptx(2);
    pt3x = ptx(3);
    pt4x = ptx(4);
        
    %find intersection point using triangulation
    xvals = cos( [pt1y,pt3y,pt2y,pt4y] .*pi/180) .* ...
        cos( [pt1x,pt3x,pt2x,pt4x] .*pi/180);
    yvals = sin( [pt1y,pt3y,pt2y,pt4y] .*pi/180) .* ...
        cos( [pt1x,pt3x,pt2x,pt4x] .*pi/180);
    zvals = sin( [pt1x,pt3x,pt2x,pt4x] .*pi/180);
    A = [xvals(1),yvals(1),zvals(1)];
    B  = [xvals(2),yvals(2),zvals(2)];
    C  = [xvals(3),yvals(3),zvals(3)];
    D =  [xvals(4),yvals(4),zvals(4)];
    N = cross(B-A,C-A);
    N = N./norm(N);
    latc1 = -asin(N(3)).*180/pi;
    lonc1= atan2(N(2),N(1))*180/pi-180;
    N2 = cross(B-A,D-A);
    N2 = N2./norm(N2);
    latc2 = -asin(N2(3))*180/pi;
    lonc2 = atan2(N2(2),N2(1))*180/pi-180;
    N3 = cross(C-B,D-B);
    N3 = N3./norm(N3);
    latc3 = -asin(N3(3)).*180/pi;
    lonc3= atan2(N3(2),N3(1))*180/pi-180;
    N4 = cross(D-C,A-C);
    N4 = N4./norm(N4);
    latc4 = -asin(N4(3))*180/pi;
    lonc4 = atan2(N4(2),N4(1))*180/pi-180;
    lat0 = mean([latc1,latc2,latc3,latc4]);
    lon0 = mean([lonc1,lonc2,lonc3,lonc4]);
    if lon0 < (-180); lon0 = lon0 + 360; end;
    if lon0 > ( 180); lon0 = lon0 - 360; end;
    lat0s(i) = lat0;
    lon0s(i) = lon0;
    latcens(i) = info.act_centroid_lat(ind(i));
    loncens(i) = -info.act_centroid_lon(ind(i));
    if loncens(i) < (-180); loncens(i) = loncens(i) + 360; end;
    if loncens(i) > ( 180); loncens(i) = loncens(i) - 360; end;
    
    %derive semi-major and semi-minor axis for geographic ellipse
    a = sqrt( (pt1x-pt2x).^2 + (pt1y-pt2y).^2 );
    b = sqrt( (pt3x-pt4x).^2 + (pt3y-pt4y).^2 );
    c = sqrt( (pt1x-pt3x).^2 + (pt1y-pt3y).^2 );
    d = sqrt( (pt1x-pt4x).^2 + (pt1y-pt4y).^2 );
    e = sqrt( (pt2x-pt3x).^2 + (pt2y-pt3y).^2 );
    f = sqrt( (pt2x-pt4x).^2 + (pt2y-pt4y).^2 );
    [major,mind] = max([a,b,c,d,e,f]);
    switch mind
        case 1
            minor = b;
            az2 = atan( (pt1y-pt2y)./(pt1x-pt2x) );
        case 2
            minor = a;
            az2 = atan( (pt4y-pt3y)./(pt4x-pt3x) );
        case 3
            minor = f;
            az2 = atan( (pt1y-pt3y)./(pt1x-pt3x) );
        case 4
            minor = e;
            az2 = atan( (pt4y-pt1y)./(pt4x-pt1x) );
        case 5
            minor = d;
            az2 = atan( (pt2y-pt3y)./(pt2x-pt3x) );
        case 6;
            minor = c;
            az2 = atan( (pt4y-pt2y)./(pt4x-pt2x) );
    end;       
    ecc = sqrt( 1 - (minor./major).^2);

    if major < small_thresh;
        yes_small = 1;
    else;
        yes_small = 0;
    end;
    
    %%find angular offsets
    %         %normalize semimajor axis
    %         n1 =[ pt1x - pt2x, pt1y-pt2y ]./sqrt( (pt1x-pt2x)^2+(pt1y-pt2y)^2);
    %         off1 = acos( n1(2) ) * 180/pi;
    %         %normalize semiminor axis
    %         n2 =[ pt3x - pt4x, pt3y-pt4y ]./sqrt( (pt3x-pt4x)^2+(pt3y-pt4y)^2);
    %         off2 = acos( n2(2) ) * 180/pi;
    
    %get the reported azimuthal look direction
    az = (-1).*(info.act_azimuth_angle(ind(i)));
    
    %create the ellipse    
    t = linspace(-pi,pi,100);
    lat = lat0 + major./2.*cos(t).*cos(az2) - ...
                 minor./2.*sin(t).*sin(az2);             
    lon = lon0 + major./2.*cos(t).*sin(az2) + ...
                 minor./2.*sin(t).*cos(az2);
                 
%     %create the ellipse
%     [lat,lon] = ellipse1(lat0,lon0,[major./2,ecc],az);

    %convert the longitude values into -180 to 180 
    lon(lon > 180) = lon(lon > 180) - 360;
    lon0s(lon0s > 180) = lon0s(lon0s > 180) - 360;
    loncens(loncens > 180) = loncens(loncens > 180) - 360;
    
    %if we pass through 180 in the middle fix it
    if ( max(lon) > 170 ) & ( min(lon) < -170 );
        lon(lon < 0) = lon(lon < 0) + 360;
    end;
    

%     %fix the sign conventions
%     loncens = -loncens;
%     lat = -lat;
%     lon = (-lon)+180;
%     lon0s(i) = -(lon0s(i))+180;
%     lat0s(i) = -lat0s(i);
       
    %populate the structure file
    struct(i).BoundingBox = ...
        [min(lat) max(lat); min(lon) max(lon)];
    struct(i).Y = lat;
    struct(i).X = lon;
    struct(i).Geometry = 'Polygon';
    struct(i).Id = i;
    struct(i).Energy = info.total_echo_energy(ind(i))*1e16;
    struct(i).Energy_Units = 'Joules * 1e16';
    struct(i).Sigma0_Cor = info.sigma0_corrected(ind(i));
    struct(i).Sigma0_UnCor = info.sigma0_uncorrected(ind(i));
    struct(i).Num_Pulse = info.num_pulses_received(ind(i));
    struct(i).Gain = info.system_gain(ind(i));
    struct(i).Antenna = info.antenna_temp(ind(i));
    struct(i).Antenna_Std = info.ant_temp_std(ind(i));
    struct(i).Receiver = info.receiver_temp(ind(i));
    struct(i).Pass_Pol_Angle = info.pass_pol_angle(ind(i));
    struct(i).Pass_Azimuth_Angle = info.pass_azimuth_angle(ind(i));
    struct(i).Pass_Emission_Angle = info.pass_emission_angle(ind(i));
    struct(i).Act_Pol_Angle = info.act_pol_angle(ind(i));
    struct(i).Act_Incidence_Angle = info.act_indcidence_angle(ind(i));
    struct(i).Act_Azimuth_Angle = info.act_azimuth_angle(ind(i));
    struct(i).X_Factor = info.x_factor(ind(i));
    struct(i).Radar_Mode = info.radar_mode(ind(i));
    struct2(i).X = loncens(i);
    struct2(i).Y = latcens(i);
    struct2(i).Geometry = 'Point';
    struct2(i).Id = i;
    struct2(i).Energy = info.total_echo_energy(ind(i))*1e16;
    struct2(i).Energy_Units = 'Joules * 1e16';
    struct2(i).Sigma0_Cor = info.sigma0_corrected(ind(i));
    struct2(i).Sigma0_UnCor = info.sigma0_uncorrected(ind(i));
    struct2(i).Num_Pulse = info.num_pulses_received(ind(i));
    struct2(i).Gain = info.system_gain(ind(i));   
    struct2(i).Antenna = info.antenna_temp(ind(i));
    struct2(i).Antenna_Std = info.ant_temp_std(ind(i));
    struct2(i).Receiver = info.receiver_temp(ind(i));
    struct2(i).Pass_Pol_Angle = info.pass_pol_angle(ind(i));
    struct2(i).Pass_Azimuth_Angle = info.pass_azimuth_angle(ind(i));
    struct2(i).Pass_Emission_Angle = info.pass_emission_angle(ind(i));
    struct2(i).Act_Pol_Angle = info.act_pol_angle(ind(i));
    struct2(i).Act_Incidence_Angle = info.act_indcidence_angle(ind(i));
    struct2(i).Act_Azimuth_Angle = info.act_azimuth_angle(ind(i));
    struct2(i).X_Factor = info.x_factor(ind(i));
    struct2(i).Radar_Mode = info.radar_mode(ind(i));
    struct3(i).X = lon0s(i);
    struct3(i).Y = lat0s(i);
    struct3(i).Geometry = 'Point';
    struct3(i).Id = i;
    struct3(i).Energy = info.total_echo_energy(ind(i))*1e16;
    struct3(i).Energy_Units = 'Joules * 1e16';
    struct3(i).Sigma0_Cor = info.sigma0_corrected(ind(i));
    struct3(i).Sigma0_UnCor = info.sigma0_uncorrected(ind(i));
    struct3(i).Num_Pulse = info.num_pulses_received(ind(i));
    struct3(i).Gain = info.system_gain(ind(i));     
    struct3(i).Antenna = info.antenna_temp(ind(i));
    struct3(i).Antenna_Std = info.ant_temp_std(ind(i));
    struct3(i).Receiver = info.receiver_temp(ind(i));
    struct3(i).Pass_Pol_Angle = info.pass_pol_angle(ind(i));
    struct3(i).Pass_Azimuth_Angle = info.pass_azimuth_angle(ind(i));
    struct3(i).X_Factor = info.x_factor(ind(i));
    struct3(i).Pass_Emission_Angle = info.pass_emission_angle(ind(i));
    struct3(i).Act_Pol_Angle = info.act_pol_angle(ind(i));
    struct3(i).Act_Incidence_Angle = info.act_indcidence_angle(ind(i));
    struct3(i).Act_Azimuth_Angle = info.act_azimuth_angle(ind(i));
    struct3(i).Radar_Mode = info.radar_mode(ind(i));
    if yes_small;
        struct4(scnt).BoundingBox = ...
            [min(lat) max(lat); min(lon) max(lon)];
        struct4(scnt).Y = lat;
        struct4(scnt).X = lon;
        struct4(scnt).Geometry = 'Polygon';
        struct4(scnt).Id = i;
        struct4(scnt).Gain = info.system_gain(ind(i));        
        struct4(scnt).Energy = info.total_echo_energy(ind(i))*1e16;
        struct4(scnt).Energy_Units = 'Joules * 1e16';
        struct4(scnt).Sigma0_Cor = info.sigma0_corrected(ind(i));
        struct4(scnt).Sigma0_UnCor = info.sigma0_uncorrected(ind(i));
        struct4(scnt).Num_Pulse = info.num_pulses_received(ind(i));
        struct4(scnt).Gain = info.system_gain(ind(i));
        struct4(scnt).Antenna = info.antenna_temp(ind(i));
        struct4(scnt).Antenna_Std = info.ant_temp_std(ind(i));
        struct4(scnt).Receiver = info.receiver_temp(ind(i));
        struct4(scnt).Pass_Pol_Angle = info.pass_pol_angle(ind(i));
        struct4(scnt).Pass_Azimuth_Angle = info.pass_azimuth_angle(ind(i));
        struct4(scnt).X_Factor = info.x_factor(ind(i));
        struct4(scnt).Pass_Emission_Angle = info.pass_emission_angle(ind(i));
        struct4(scnt).Act_Pol_Angle = info.act_pol_angle(ind(i));
        struct4(scnt).Act_Incidence_Angle = info.act_indcidence_angle(ind(i));
        struct4(scnt).Act_Azimuth_Angle = info.act_azimuth_angle(ind(i));
        struct4(scnt).Radar_Mode = info.radar_mode(ind(i));
        scnt = scnt+1;
    end;
end;

%write out the shape file
shapewrite(struct,file2_out);
shapewrite(struct2,file2_out2);
shapewrite(struct3,file2_out3);
shapewrite(struct4,file2_out4);

%create projection files for the new shape files
%open and write a projection file for titan to go along with it
out_file = file2_out;
pos = strfind( out_file, '.shp' );
prj_file = [out_file(1:pos) 'prj'];
fid = fopen(prj_file,'w');
string = ['GEOGCS["GCS_Titan_2000",DATUM["D_Titan_2000",'...
    'SPHEROID["Titan_2000_IAU_IAG",2575000.0,0.0]],PRIMEM["Reference_Meridian",0.0],',...
    'UNIT["Degree",0.0174532925199433]]' '\n'];
fprintf(fid,string);
fclose(fid);
out_file = file2_out2;
pos = strfind( out_file, '.shp' );
prj_file = [out_file(1:pos) 'prj'];
fid = fopen(prj_file,'w');
string = ['GEOGCS["GCS_Titan_2000",DATUM["D_Titan_2000",'...
    'SPHEROID["Titan_2000_IAU_IAG",2575000.0,0.0]],PRIMEM["Reference_Meridian",0.0],',...
    'UNIT["Degree",0.0174532925199433]]' '\n'];
fprintf(fid,string);
fclose(fid);
out_file = file2_out3;
pos = strfind( out_file, '.shp' );
prj_file = [out_file(1:pos) 'prj'];
fid = fopen(prj_file,'w');
string = ['GEOGCS["GCS_Titan_2000",DATUM["D_Titan_2000",'...
    'SPHEROID["Titan_2000_IAU_IAG",2575000.0,0.0]],PRIMEM["Reference_Meridian",0.0],',...
    'UNIT["Degree",0.0174532925199433]]' '\n'];
fprintf(fid,string);
fclose(fid);
out_file = file2_out4;
pos = strfind( out_file, '.shp' );
prj_file = [out_file(1:pos) 'prj'];
fid = fopen(prj_file,'w');
string = ['GEOGCS["GCS_Titan_2000",DATUM["D_Titan_2000",'...
    'SPHEROID["Titan_2000_IAU_IAG",2575000.0,0.0]],PRIMEM["Reference_Meridian",0.0],',...
    'UNIT["Degree",0.0174532925199433]]' '\n'];
fprintf(fid,string);
fclose(fid);

disp(['Done with Scatterometry For ' filename]);
end;

disp(['Done with ' filename]);