function out = GetSRTMData_new(p1, pathin)
% function out = GetSRTMData(p1, [p2])
%
% General:
% The function retrieves SRTM-Data from the according SRTM-files (3 arcsecond resolution only !),
% which need to be supplied by the user. The path to the unzipped SRTM-files
% (.hgt-format) needs to be specified in the variable 'path'.
% SRTM-files can be downloaded for free -> ftp://e0dps01u.ecs.nasa.gov/srtm/
% Missing data can be filled with the free tool SRTM-Fill  -> http://www.3dnature.com/srtmfill.html
%
% Input Argument(s):
% p1/2 need to be specified as decimal degree coordinates,
% negative values for western longitudes resp. southern latitudes.
% p1  			-  lon/lat point or lon/lat point-list
% p2(optional)  -  lon/lat point or lon/lat point-list (list of same size as p1)
%
% Output Argument:
% out - structure containing Point(s)/profiles with [lon / lat / height].
% ...
% The following height data is retrived from the SRTM-files:
% p1 point or list / p2 empty: -> For specified point(s) 							-> Example1
% p1 point / p2 point: -> For points on line p1 <-> p2 								-> Example2
% p1 list / p2 list: -> For points on line p1(1) <-> p2(1) ... p1(n) <-> p2(n)		-> Example3
%
% Examples:
% 1: out = GetSRTMData([100.5 40.1234]);
% 2: out = GetSRTMData([100.5 40.1234],[100.6 41.4321]);
% 3: out = GetSRTMData([100 40; 101 41];[100.1 40.8; 101.1 41.8]);
%
% Version 0.9:  Just finished it.
% Version 0.91: Added switch for files with southern latitudes / western longitudes.
% Version 1.0: Corrected bug mentioned by John Tanner for lat / lon <10° (3.12.2004).
%
% Author
% Sebastian Hölz, TU Berlin, Germany
% hoelz@geophysik.tu-berlin.de
% If you find any bugs, let me know ...
p2 = [];
if nargin ==1; p2=[]; end

% 1. Checking for right dimension (nx2)
if size(p1,2) ~= 2; p1=p1'; end
if size(p2,2) ~= 2; p2=p2'; end
if nargin == 3 && size(p1) ~= size(p2)
    m = msgbox('For height profiles the point lists need to have the same size');
    uiwait(m);
    out = [];
    return
end

% 1. Loading SRTM-files
[SRTM, lon, lat, abort, out] = LoadSRTM(p1,p2,pathin);
if abort
    out = [];
    return;
elseif and(~abort,~isempty(out))
    return;
end

% 2. Generating the profile indizes
if isempty(p2)
    out = {GetIndizes(lon, lat, p1, p2)};
else
    for i = 1:size(p1,1)
        out = GetIndizes(lon, lat, p1(i,:), p2(i,:));
    end
end

% 3. Extracting Altitude data
for i = 1:length(out)
    ind = (out{i}(:,1)-1)*size(SRTM,1)+out{i}(:,2);
    out = SRTM(ind);
    %out{i}(:,1)=(out{i}(:,1)-1)/1200+lon(1);
    %out{i}(:,2)=(out{i}(:,2)-1)/1200+lat(1);
end

% if length(out)==1; out=out{1}; end

% -----------------------------------------
function [SRTM, lon, lat, abort, out] = LoadSRTM(p1,p2,pathin)

global missing_filenames

out = [];
abort = 0; % does the program abort in this subfunction (0 = no, 1 = yes )
SRTM=[];
SRTM_tmp=[];
regions = {'Eurasia','North_America','South_America','Afrika','Islands'}; % possible continent-url-extensions
tmp = [p1; p2];

lon(1) = floor(min(tmp(:,1)));
lon(2) = ceil(max(tmp(:,1)));
lat(1) = floor(min(tmp(:,2)));
lat(2) = ceil(max(tmp(:,2)));
if lon(2)==lon(1); lon(2)=lon(2)+1; end
if lat(2)==lat(1); lat(2)=lat(2)+1; end

for x = lon(1):lon(2)-1
    for y = lat(1):lat(2)-1
        if x<=-100; LON = 'W';
        elseif x>-100 && x<=-10; LON = 'W0';
        elseif x>-10 && x<0; LON = 'W00';
        elseif x>=0 && x<10; LON = 'E00';
        elseif x>=10 && x<100; LON = 'E0'; %NEW: before: elseif x>10 & x<100; LON = 'E0'--> mentioned as a Bug in FEX (post january 2011)
                                           % --> thats right, because I had this case as well (program didn't find available SRTM-file because of wrong naming)
        else LON ='E';
        end
        if y<=-10; LAT = 'S';
        elseif y>-10 && y<0; LAT = 'S0';
        elseif y>=0 && y<10; LAT = 'N0';
        else LAT = 'N';
        end
        SRTMname = [LAT num2str(abs(y)) LON num2str(abs(x))];
        file=[pathin 'unzipped' filesep SRTMname '.hgt'];
        fid=fopen(file,'r','b');
        while fid == -1
            for k = 1:length(regions)
                [junk,status] = urlwrite(['http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/' regions{k} '/' SRTMname '.hgt.zip'], ...
                    [pathin SRTMname '.hgt.zip']);
                if status
                    break;
                elseif and(~status,k==length(regions)) % no regions left to try in the URL
                    if ~any(strcmp(missing_filenames,SRTMname)) % Is the current SRTM-file NOT in the missing_filenames variable included?
                        % --> then ask the user if he would like to try to download the elevation data from http://ws.geonames.org
                        m = questdlg(sprintf('%s\n%s',['Required SRTM-file (' SRTMname ') is not online available.'],...
                            'Try to get the elevation data from http://ws.geonames.org for this file?'),'?','yes','no','yes');
                        pause(0.1);
                        if isempty(m) || strcmp(m,'no')
                            abort = 1;
                            return;
                        end
                        % User pressed 'yes' OR the current coordinates still belong to an already as "unavailabe" known SRTM-file (saved in Cell "missing_filenames")
                        % --> the second point avoids to ask the user every run if he would like to download elevation data
                        % from http://ws.geonames.org for this SRTM-files
                        
                        if isempty(missing_filenames)
                            missing_filenames{1} = SRTMname;
                        else
                            missing_filenames{end+1} = SRTMname;
                        end
                    end
                    loopnum = 0;
                    while 1
                        loopnum = loopnum+1;
                        elevation = str2num(urlread(['http://ws.geonames.org/srtm3?lat=' num2str(tmp(2)) '&lng=' num2str(tmp(1))]));
                        if and(isnumeric(elevation),~isempty(elevation)) % Is elevation a numeric value-->otherwise the request failed
                            out = elevation; 
                            return;
                        elseif loopnum == 20 % Abort the request if the server did not response after x-attempts (e.g. 20) with a valid value
                            m = msgbox('Server overloaded: Could''t get elevatin data','error');
                            uiwait(m);
                            abort = 1;
                            return;
                            
                        end
                    end
                end
            end
            
            unzip([pathin SRTMname '.hgt.zip'],[pathin 'unzipped']);
            fid = fopen(file,'r','b');
        end
        
        if isempty(SRTM_tmp)
            SRTM_tmp = rot90(fread(fid,[1201 1201],'*int16')); % NEW: changed from uint16 --> mentioned in FEX that there are also towns under the sea level 
        else
            SRTM_tmp = [SRTM_tmp(1:end-1,:); rot90(fread(fid,[1201 1201],'*int16'))]; %NEW: same here
        end
        fclose(fid);
    end
    if isempty(SRTM)
        SRTM = SRTM_tmp;
    else
        SRTM = [SRTM(:,1:end-1) SRTM_tmp];
    end
    SRTM_tmp = [];
end

% ------------------------------
function out = GetIndizes(lon, lat, p1, p2)

if isempty(p2)                                  % Only single point or point list needed
    out = zeros(size(p1,1),3);
    out(:,1)=round((p1(:,1)-lon(1))*1200)+1;
    out(:,2)=round((p1(:,2)-lat(1))*1200)+1;
else                                            % Profile between two points needed
    indx = [round((p1(1)-lon(1))*1200)+1 round((p2(1)-lon(1))*1200)+1];
    indy = [round((p1(2)-lat(1))*1200)+1 round((p2(2)-lat(1))*1200)+1];
    diffx = diff(indx);
    diffy = diff(indy);
    len = max(abs([diffx diffy]))+1;
    
    out = zeros(len,2);
    if indx(1)==indx(2)
        out(:,1)=indx(1);
    else
        out(:,1)=round(indx(1):diffx/(len-1):indx(2))';
    end
    if indy(1)==indy(2)
        out(:,2)=indy(1);
    else
        out(:,2)=round(indy(1):diffy/(len-1):indy(2))';
    end
end