function [pyy, f] = psdexamplescript2(strFile, fs, Nfft, No, w)

% % EXAMPLE
% strURL = 'https://gipoc.grc.nasa.gov/pims/pub/pad/year2022/month05/day01/sams2_accel_121f08/2022_05_01_01_01_03.799+2022_05_01_01_11_03.804.121f08';
% strFile = '/tmp/2022_05_01_01_01_03.799+2022_05_01_01_11_03.804.121f08'; % download file above and save it to local machine tmp directory
% urlwrite(strURL, strFile); % you may not have this function, but can download by whatever means and save the file locally (in say tmp directory)
% fs = 500; % SampleRate in sa/sec; see [strURL '.header']
% Nfft = 16384;
% No = Nfft/2;
% w = Nfft; % window; be careful in MATLAB version with how this parameter is handled
%           % for PIMS routines, w = Nfft is typical usage and implies hanning(Nfft)
% [pyy, f] = psdexamplescript2(strFile, fs, Nfft, No, w);
%
% NOTES -- see https://gipoc.grc.nasa.gov/pims/plots/batch/year2022/month05/day01/2022_05_01_00_00_00.000_121f08_spgs_roadmaps500.pdf
% - a PSD just after GMT hour zero would show no signs of strong narrowband spectral peaks around 60 Hz
% - a PSD just after GMT hour one would shows strong signs of narrowband disturbance (spectral peaks) around/above 60 Hz

% Load data from binary file, and ...
% note from [strURL '.header'] file that these data were archived in SSA coordinates
fid = XXXXX(strFile, 'r', 'l');
data = XXXXX(fid, [4 inf], 'float')'; % for MAMS, need to change 4 to number of columns in binary PAD file
fclose(fid);
y = data(:, 3); % data is [ t x y z ], so y-axis is 3rd column

% Compute PSD for SSA Y-axis (call it pyy)
[pyy, f] = pwelch(y, w, No, Nfft, fs);

% Crude verify via Parseval
rms_td = sqrt(mean(y.^2));      % RMS computed from signal in time domain
rms_fd = sqrt(trapz(f,pyy));    % RMS computed from signal in freq domain
fprintf('A crude PSD verification via Parseval''s Theorem:\nRMS from signal in time domain = %e\nRMS from signal in freq domain = %e\n', rms_td, rms_fd)

% The beauty of Parseval (integrate over frequency band(s) of interest)
% for this data set shows that a good deal of RMS comes from the spectral
% peaks around 60 Hz (between GMT hour one and two)
%ugrms_bands = parseval(f, pyy, 200, [0 55; 55 65; 65 200]) / 1e-6;  % 3 frequency bands: compute RMS in each band via Parseval
