function specData = getDesignSpectrumASCE7_22_LatLon(latitude, longitude, siteClass, riskCategory)
% HARDENED ASCE 7-22 USGS QUERY (LAT/LON VERSION)
%
% Features:
%   - curl (IPv4 only) → eliminates TLS/IPv6 failures
%   - retry with exponential backoff
%   - response validation
%   - caching
%   - fallback if offline
%
% Author: Tom Irvine (fully hardened)

%% -----------------------------
% Input handling
%% -----------------------------
siteClass = upper(strtrim(char(siteClass)));
[riskRoman,~] = convertRiskCategory(riskCategory);

lat = latitude;
lon = longitude;

%% -----------------------------
% Cache setup
%% -----------------------------
cacheDir = fullfile(tempdir,'usgs_cache');
if ~exist(cacheDir,'dir')
    mkdir(cacheDir);
end

cacheFile = fullfile(cacheDir, ...
    sprintf('lat%.4f_lon%.4f_%s_%s.json',lat,lon,siteClass,riskRoman));

%% -----------------------------
% Try cache first
%% -----------------------------
d22 = [];

if exist(cacheFile,'file')
    try
        txt = fileread(cacheFile);
        resp = jsondecode(txt);
        d22 = resp.response.data;
        disp('USGS data loaded from cache');
    catch
        d22 = [];
    end
end

%% -----------------------------
% Query USGS if needed
%% -----------------------------
if isempty(d22)

    base22 = 'https://earthquake.usgs.gov/ws/designmaps/asce7-22.json';
    titleParam = 'LatLonQuery';

    url22 = sprintf(['%s?latitude=%.6f&longitude=%.6f', ...
                     '&riskCategory=%s&siteClass=%s&title=%s'], ...
                     base22,lat,lon,riskRoman,siteClass,titleParam);

    resp = [];

    for attempt = 1:5

        fprintf('USGS query attempt %d...\n',attempt);

        try
            fname = [tempname '.json'];

            cmd = sprintf(['curl -4 -L --retry 2 --retry-delay 2 ', ...
                           '--connect-timeout 10 --max-time 30 ', ...
                           '"%s" -o "%s"'], url22, fname);

            [status,~] = system(cmd);

            if status==0 && exist(fname,'file')
                txt = fileread(fname);
                tmp = jsondecode(txt);

                if isfield(tmp,'response') && isfield(tmp.response,'data')
                    resp = tmp;
                    fprintf('✔ curl success\n');
                    break
                end
            end
        catch
        end

        pause(2^attempt);  % exponential backoff
    end

    %% -----------------------------
    % fallback if all fail
    %% -----------------------------
    if isempty(resp)

        warning('USGS unavailable → using fallback values');

        Ss = 1.5;
        S1 = 0.6;
        TL = 8;
        periods = [];
        ordinates = [];

    else

        % Save cache
        try
            fid = fopen(cacheFile,'w');
            fwrite(fid,jsonencode(resp),'char');
            fclose(fid);
        end

        d22 = resp.response.data;

        Ss = getFieldSafe(d22,'ss');
        S1 = getFieldSafe(d22,'s1');
        TL = getFieldSafe(d22,'tl');

        if isfield(d22,'multiPeriodDesignSpectrum')
            periods   = d22.multiPeriodDesignSpectrum.periods;
            ordinates = d22.multiPeriodDesignSpectrum.ordinates;
        else
            periods   = [];
            ordinates = [];
        end
    end

else

    Ss = getFieldSafe(d22,'ss');
    S1 = getFieldSafe(d22,'s1');
    TL = getFieldSafe(d22,'tl');

    if isfield(d22,'multiPeriodDesignSpectrum')
        periods   = d22.multiPeriodDesignSpectrum.periods;
        ordinates = d22.multiPeriodDesignSpectrum.ordinates;
    else
        periods   = [];
        ordinates = [];
    end
end

%% -----------------------------
% Compute Fa, Fv
%% -----------------------------
Fa = computeFa(Ss,siteClass);
Fv = computeFv(S1,siteClass);

SMS = Fa*Ss;
SM1 = Fv*S1;

SDS = (2/3)*SMS;
SD1 = (2/3)*SM1;

%% -----------------------------
% Output
%% -----------------------------
specData = struct( ...
    'siteClass',siteClass, ...
    'riskCategory',riskRoman, ...
    'latitude',lat, ...
    'longitude',lon, ...
    'Ss',Ss,'S1',S1, ...
    'Fa',Fa,'Fv',Fv, ...
    'SMS',SMS,'SM1',SM1, ...
    'SDS',SDS,'SD1',SD1, ...
    'TL',TL, ...
    'periods',periods, ...
    'ordinates',ordinates);

end

%% =============================================================
% Helper Functions
%% =============================================================

function [riskRoman, riskNumeric] = convertRiskCategory(riskCategory)

if isnumeric(riskCategory)
    riskNumeric = riskCategory;
else
    s = upper(strtrim(char(riskCategory)));
    switch s
        case 'I',   riskNumeric = 1;
        case 'II',  riskNumeric = 2;
        case 'III', riskNumeric = 3;
        case 'IV',  riskNumeric = 4;
        otherwise
            riskNumeric = str2double(s);
    end
end

if ~ismember(riskNumeric,[1 2 3 4])
    error('Risk category must be I–IV or 1–4.');
end

switch riskNumeric
    case 1, riskRoman = 'I';
    case 2, riskRoman = 'II';
    case 3, riskRoman = 'III';
    case 4, riskRoman = 'IV';
end
end

function val = getFieldSafe(s,fld)
if isstruct(s) && isfield(s,fld)
    val = s.(fld);
else
    val = NaN;
end
end

function Fa = computeFa(Ss, siteClass)

switch siteClass
    case 'A', Fa = 0.8;
    case 'B', Fa = 1.0;
    case 'C'
        Fa = interp1([0.25 0.50 0.75 1.00], ...
                     [1.2 1.1 1.0 1.0], ...
                     min(max(Ss,0.25),1.0),'linear');
    case 'D'
        Fa = interp1([0.25 0.50 0.75 1.00], ...
                     [1.6 1.4 1.2 1.1], ...
                     min(max(Ss,0.25),1.0),'linear');
    case 'E'
        Fa = interp1([0.25 0.50 0.75 1.00], ...
                     [2.5 1.7 1.2 0.9], ...
                     min(max(Ss,0.25),1.0),'linear');
    otherwise
        Fa = 1.0;
end
end

function Fv = computeFv(S1, siteClass)

switch siteClass
    case 'A', Fv = 0.8;
    case 'B', Fv = 1.0;
    case 'C'
        Fv = interp1([0.1 0.2 0.3 0.4], ...
                     [1.7 1.6 1.5 1.4], ...
                     min(max(S1,0.1),0.4),'linear');
    case 'D'
        Fv = interp1([0.1 0.2 0.3 0.4], ...
                     [2.4 2.0 1.8 1.6], ...
                     min(max(S1,0.1),0.4),'linear');
    case 'E'
        Fv = interp1([0.1 0.2 0.3 0.4], ...
                     [3.5 3.2 2.8 2.4], ...
                     min(max(S1,0.1),0.4),'linear');
    otherwise
        Fv = 1.0;
end
end