 function [noise_ASD, response_FR, response_ZPK] = SEI_sensor_noise(sensor_name, freq,verbose)
% SEI_sensor_noise  Noise estimates for HAM-ISI, BSC-ISI, and HEPI SEI sensors
%  version 4
%
%  [noise, vargargout] = SEI_sensor_noise('sensor_name', freq, varargin)
%  noise is the ASD of the self noise of the sensor, in m/rtHz
%    this includes the noise of the default LIGO preamp, if any
%    see T0900450 for more details.
%  freq is the frequency vector (Hz)
%  sensor_name is a string defining the amplifier in question
%   so far 
%          'ADE_1mm', 'ADE_p25mm', 'Kaman_IPSmeas',     % Displacement Sensors
%          'Fine_Hugh', 
%          'L4C', 'GS13meas', 'GS13calc',               % 1 Hz Geophones
%          'T240meas', 'T240spec', 'Guralp_40T_Spec',   % Low-freq Seismometers 
%          'T120H' 'T120C' 'T360'
%          'BRS', 'cBRS', 'CRS'                         % Rotation Sensors (rad/rtHz)
%          'PCB_356M132_spec', 'Wilcoxon_797L_spec'     % Accelerometers
%          'PICO', 'SmarAct', 'HoQI', 'SPI_L',          % Interferometric Sensors
%          'SPI_PY'                                     % Optical Levers (rad/rtHz)
%   have been included as well as 
%          'GS_16bit_ADC'    which returns in Vrms/rtHz, not m/rtHz.
%   in 2021 BTL added
%          'projections', which returns a data structure to project the
%          noise of 1 sensor to various cartesian base. Multiply the noise
%          of 1 sensor by the correct element in this struct to get the
%          noise for the CART DOF. see T2100336
%
% (the code is not case sensitive)
% This can also be called with an optional second output argument
% which returns the frequency response of the sensitivity of the sensor
% in volts/ meter, e.g.
%
% [noise, response_FR] = SEI_sensor_noise('GS13meas',freq)
% returns the ASD of displacement noise (noise_ASD) for the GeoTech GS13
% 'GS13' as well as the frequency response of the instrument (response_FR)
%
% noise = SEI_sensor_noise('L4C',freq,true)
% returns the ASD of displacement noise (noise_ASD) only, but also spits
% out several comments about the noise when the third boolean argument
% (verbose) is set to true.
%
% 2009-09-24 Authored Brian Lantz, 
% 2010-03-26 ADC noise added from Jay Heefner
% 2011-10-18 BTL update the T240 response
% 2014-02-25 Added factor of two to noise of 1mm CPS ( noise_ASD    = 2*10.^lognoise;)  
% 2014-11-05 Added Guralp  
% 2017-02-03 Added PCB and Wilcoxon Accelerometers, updated help text and
%            sensor input options 
% 9/2018 RKM
% added T360 from page 3 of https://www.nanometrics.ca/sites/default/files/2018-04/trillium360_brochure_.pdf
% added PICO PICOscale interfermoter from
% http://www.smaract.com/wp-content/uploads/2017/01/PS-SS00007_ResolutionAndStability-1.pdf
% and other communication 
 % 4/2019 RKM
 % updated T360 and T240 to use data from Nanometrics (current T360)
 % added T120H (Horizon) and T120C (compact) with data from nanometrics
 % 3/16/20
 %added a CPS model called Fine_Hugh which matchs the data in SEI 1580,
 % looks like this http://prntscr.com/ruyz5e
 % 8/3/2021 - BTL added the noise projections. If you call as 
 %   SEI_sensor_noise('projections') you will get back the data structure
 %   for the projections from individual sensor noise to the cartesian basis.
 % 6/30/2023 -- JSK added cBRS, CRS, HoQI, SPI_L, and SPI_P.
 
  if nargin < 3  %default to not verbose
      verbose = 0;
  end
  
  if verbose
      disp(' ')
      disp('you are using draft version 5. call Brian Lantz for the real one!')
      disp(' ')
  end

if strncmpi(sensor_name,'ADE_1mm',7)
    % see Andy Stein SEI log enty 1311, nov 4 overnight plot
    freq_data   = [.001,  .002,  .01,    .1,    .7,    100];
    noise_data  = [ 5e-8, 9e-9, 2.5e-9, 5e-10, 2e-10, 2e-10];
    
    lognoise     =  interp1(log10(freq_data),log10(noise_data),log10(freq));
    noise_ASD    = 2*10.^lognoise;
    response_FR  = 10/1e-3 * ones(size(freq)); % 10 V per mm
    response_ZPK = tf(10/1e-3);


elseif strncmpi(sensor_name,'Fine_Hugh',9)
    % aPLUS order for Fine CPS (ADE_p25mm) that are better than the aLIGO
    % order; see SEI aLOG 1580
   
      A = log( 5.5e-10/2e-10)/log(0.03/0.1);
      B = log(5.5e-10) - A*log(0.03);
      N = exp(A*log(0.001) +B);
      freq_data   =   [0.001        .03,    0.1,         1,        2.5,      100 1000];
    noise_data  = [     N     5.5e-10,     2e-10       0.5e-10,     4e-11,     3e-11 3e-11;];
%     
      lognoise    = interp1(log10(freq_data),log10(noise_data),log10(freq));
      noise_ASD   = 10.^lognoise;
%     response_FR = 10/0.25e-3 * ones(size(freq)); % 10V/ 0.25 mm
%     response_ZPK = tf(10/0.25e-3); 
      response_FR =  ones(size(freq)); % 10V/ 0.25 mm
     response_ZPK = tf(1);    
     
elseif strncmpi(sensor_name,'ADE_p25mm',9)
	    freq_data   =   [.001,   .002,    .01,      .1,    .7,    100];
	    noise_data  = [1.5e-08, 2.7e-09, 7.5e-10, 1.5e-10, 6e-11, 6e-11];
        lognoise    = interp1(log10(freq_data),log10(noise_data),log10(freq));
        noise_ASD   = 10.^lognoise;
         response_FR =  ones(size(freq)); % 10V/ 0.25 mm
    response_ZPK = tf(1); 
 

elseif strncmpi(sensor_name,'Kaman_IPSmeas',5)
    % See E0900426
    freq_data   = [0.001    .01,      .05,   0.1,      0.5,  1,    10,   100 ];
    noise_data  = [1e-7   1.1e-08,   1e-9, 4.5e-10, 3e-10 3e-10, 3e-10, 3e-10];
    
    lognoise    = interp1(log10(freq_data),log10(noise_data),log10(freq));
    noise_ASD   = 10.^lognoise;
    response_FR = 390/1e-3 * ones(size(freq)); % 390 V / mm
    response_ZPK = tf(390/1e-3);
    
elseif strncmpi(sensor_name,'L4C',3)
    freq_data   =   [1E-4 .04,   .52,     .8,       1.4,     4,       10,      100,  1000];
    noise_data  = [1E3 1.0e-06, 1.0e-10, 2.3e-11, 7.0e-12, 2.3e-12, 8.5e-13, 8.0e-14, 7.0e-15];
    % Extrapolate to 1k--HughR
    lognoise    = interp1(log10(freq_data),log10(noise_data),log10(freq));
    noise_ASD   = 10.^lognoise;
    
    response_ZPK = tf(0);  % this is clearly wrong, BTL Oct 18, 2011
    response_FR  = squeeze(freqresp(response_ZPK, 2*pi*freq));
    
elseif strncmpi(sensor_name,'T240spec',5)
    % this is taken directly from the T240_noise_spec_ASD.m code
%     data = [...
%         0.001	-171.688; ...
%         0.003	-179.481; ...
%         0.007	-183.636; ...
%         0.018	-187.532; ...
%         0.044	-189.610; ...
%         0.082	-190.130; ...
%         0.226	-190.130; ...
%         0.530	-189.091; ...
%         1.000	-187.273; ...
%         2.257	-183.377; ...
%         3.634	-180.000; ...
%         6.335	-174.286; ...
%         9.803	-168.312; ...
%         100     -138];        % Last line added by JSK 2012-05-02, Makes 
                              % displacement spectra role of as 1/f^(1/2), 
                              % 'cause that's what it looked like it was 
                              % turning into at 10 Hz. Complete guess, and
                              % it doesn't matter, since we never use the
                              % T240s above 10Hz anyways; it's just to get
                              % noise coupling predictions up to 100 Hz to
                              % not be identically zero
                              
                       %now using data from nanometrics         
data = [0.0001	-146.9
0.00032	-161.9
0.001	-172.1
0.0032	-179.4
0.01	-185.7
0.032	-188.9
0.1	-190.6
0.32	-190.2
1	-187.6
3.2	-180.6
10	-168.2
32	-145.8
56	-131.2
100	-116.7];
                              
    
    Freq = data(:,1);
    ANdb    = data(:,2);
      
     Accnoisedb = interp1(log(Freq), ANdb,log(freq),'pchip');
     Accnoise_ASD = 10.^(Accnoisedb/20);
     noise_ASD = Accnoise_ASD./(2*pi*freq).^2;
    
   
    % original response, changed to resp per manual on Oct 18, 2011, BTL
    %response_FR = 1196 .* 2*pi*freq;  % V/m 
    %disp('warning, the T240 response is only valid between about 20 mHz and 40 Hz')
    % 
    % the new TF comes from pg 10 of the Trillium 240 OBS user guide
    % originally specified in RAD/SEC, with minus signs. 
    %hence the lack of the usual -2*pi*
    spec_k = 2.316e9; % from the manual
    spec_gain = 1196.5;
    
    % this is in units of V/(m/s);
    response_zpk_temp = spec_gain * zpk([0 0 -108 -161], ...
        [-0.01815 + 0.01799i, -0.01815 - 0.01799i, -173, ...
        -196 + 231i, -196 - 231i, -732 + 1415i,  -732 - 1415i], spec_k) ;
    
    % and this is in V/m
    response_ZPK = tf('s') * response_zpk_temp;
    response_FR = squeeze(freqresp(response_ZPK, 2*pi*freq));
    
 
elseif strncmpi(sensor_name,'T240meas',5)
    %eyeballed from the 090310 data sets
    if verbose
        disp('The 0.01 to 0.03 Hz data from T240meas is just 2*T240spec, and not really measured')
        disp('    (the 0.03 Hz point is measured)')
    end
    freq_data   = [1E-4 .01, .03,    .1,    0.3,    1,     3,     10   100];
    noise_data  = [0.01 2*1.4e-7, 2e-8, 1.5e-9, 2e-10, 4e-11, 5e-12, 1.5e-12 1e-13];
     
    lognoise    = interp1(log10(freq_data),log10(noise_data),log10(freq));
    noise_ASD   = 10.^lognoise;

    % original response, changed to resp per manual on Oct 18, 2011, BTL
    %response_FR = 1196 .* 2*pi*freq;  % V/m 
    %disp('warning, the T240 response is only valid between about 20 mHz and 40 Hz')
    % 
    % the new TF comes from pg 10 of the Trillium 240 OBS user guide
    % originally specified in RAD/SEC, with minus signs. 
    %hence the lack of the usual -2*pi*
    spec_k = 2.316e9; % from the manual
    spec_gain = 1196.5;
    
    % this is in units of V/(m/s);
    response_zpk_temp = spec_gain * zpk([0 0 -108 -161], ...
        [-0.01815 + 0.01799i, -0.01815 - 0.01799i, -173, ...
        -196 + 231i, -196 - 231i, -732 + 1415i,  -732 - 1415i], spec_k) ;
    
    % and this is in V/m
    response_ZPK = tf('s') * response_zpk_temp;
    response_FR = squeeze(freqresp(response_ZPK, 2*pi*freq));

elseif strncmpi(sensor_name,'GS13meas',5)
    % taken directly from GS13_noise_measured_March2007 
    % the noise floor of the GS13 measured on the Tech Demo
    %  this is based on the data from the Tech Demo 3/15/2007
    %  using 2 witnesses with good ADCs, the new readouts, and a Q of ~ 5
    %  it is a little bigger than the expected noise based on the
    %  spec sheet for the LT1012 readout.
    %
    %  BTL, Oct 5, 2008
    %
    % see log entry http://ligo.phys.lsu.edu:8080/SEI/1288

    % data from graphclicks
    
    % at low freq, it looks to scale as x10 for every x2 in freq
    % is a power law of -3.32 (expect -3.5 from 1/f and sens scaling)
    % we will use the 3.5 to be conservative
    % 2e-8 * 10^3.5 = 6.3e-5
    % also assume a 1/f falloff at high freq
    full_data = [...
        1E-4   670;...
        0.01   2e-8 * 10^3.5;...
        0.101	2.138e-8;...
        0.201	1.958e-9;...
        0.400	1.892e-10;...
        0.792	1.732e-11;...
        0.994	7.262e-12;...
        1.258	3.887e-12;...
        1.655	2.449e-12;...
        2.588	1.462e-12;...
        3.954	9.725e-13;...
        8.202	4.423e-13;...
        10.710	3.280e-13;...
        24.796	1.137e-13;...
        46.087	5.610e-14;...
        91.208	2.053e-14
        1e3      2e-15;];
    
    data_freq   = full_data(:,1);
    data_noise  = full_data(:,2);
    lognoise   = interp1(log10(data_freq),log10(data_noise),log10(freq));
    noise_ASD  = 10.^lognoise;

    response_ZPK = tf(0);
    response_FR = NaN * ones(size(freq));
    disp([' warning, the freq resp for the ' sensor_name ' is not yet defined']);
    
elseif strncmpi(sensor_name,'GS13calc',6)
    freq_data   =   [.01    .1       .5      .8        1.2     3       10      100];
    noise_data  = [1.9e-5, 6.1e-9, 1.8e-11, 3.0e-12, 1.3e-12, 5.3e-13, 1.4e-13, 1.3e-14];
    
    lognoise    = interp1(log10(freq_data),log10(noise_data),log10(freq));
    noise_ASD   = 10.^lognoise;

    response_ZPK = tf(0);
    response_FR = NaN*ones(size(freq));
    disp([' warning, the freq resp for the ' sensor_name ' is not yet defined']);

elseif strncmpi(sensor_name,'GS_16bit_ADC',3)
    % General Standards 16 bit ADC PCI66-16AI64SSA
    % http://www.generalstandards.com/specs/66_16ai64ssa_c_spec_100710.pdf
    adcGain = 2^16/40; % [ct/V]
     % 0.5 [mHz] to 100 [Hz] data from Keith Riles, pg 44&45 of G1300997.
    % 1000 [Hz] data point is a guess that it stays flat.
     freq_data           = [0.0005,  0.001,   0.01,    0.1,      1,      10,  100, 1000 , 10000];
    noise_data_ct_rtHz  = [5e-1,   3.9e-1, 9.5e-2, 3.1e-2, 1.05e-2, 6.5e-3, 6e-3, 6e-3, 6e-3]; % [ct/rtHz]
    noise_data_V_rtHz   = noise_data_ct_rtHz / adcGain; % [V/rtHz]
    
    lognoise    = interp1(log10(freq_data),log10( noise_data_V_rtHz),log10(freq));
    noise_ASD   = 10.^lognoise;

    response_ZPK = tf(65536/40);
    response_FR  = squeeze(freqresp(response_ZPK, 2*pi*freq));
    disp(['          warning, the output for ' sensor_name ' has units of Volts/rtHz (not m/rtHz)']);
    
elseif strncmpi(sensor_name,'Guralp_40T_Spec',6)
    % Spec taken from data sheet found on paper in PEM drawer at LHO
    % Spec is assumed to be the bottom of the dynamic range, and captured
    % by-eye. I've extrapolated the noise down to 0.01 [Hz] and out to 100
    % [Hz], assuming the slope at either end of the spec does not change.
    period_data         = [0.01 0.02 0.03 0.04 0.07 0.11 0.25    1   10   50   100]; % [s]
    noise_data_accel_dB = [-90  -110 -121 -128 -140 -150 -160 -168 -165 -162  -161]; % [nonsense units]
    
    freq_data  = fliplr(1./period_data);
    bin_width   = (2^(1/6) - 2^(-1/6))*freq_data; % [Hz]. 1/3 octave binwidth, around central frequency of freq_data
    accel_mps2_pkpk = (10.^(fliplr(noise_data_accel_dB)/20)); % [m/s^2, pkpk] "equal to 6 times the rms acceleration for a 1/3rd octave bandwidth"
    accel_mps2_rms_prtHz = accel_mps2_pkpk ./ (6*sqrt(bin_width)); % [(m/s^2)/rtHz]
    noise_data = accel_mps2_rms_prtHz ./ ((2*pi*freq_data).^2); % [m/rtHz]
    
    lognoise   = interp1(log10(freq_data),log10(noise_data),log10(freq));
    noise_ASD  = 10.^(lognoise);
    
    response_ZPK = tf(0);
    response_FR  = squeeze(freqresp(response_ZPK, 2*pi*freq));
    disp([' warning, the freq resp for the ' sensor_name ' is not yet defined']);
    
elseif strncmpi(sensor_name,'PCB_356M132_spec',6)
    % http://www.pcb.com/Linked_Documents/Aerospace/AD_Low_Outgassing_Lowres.pdf
    % PCB Piezotronics model 356M132
    % 500 mV / g = 0.05102 V / (m/s^2)
    % 0.4 Hz to 4 kHz
    
    freq_data = 10.^(0:4);
    accel_ug_prtHz = [18 6 2 0.7 0.2];
    disp_m_prtHz   = 9.8e-6 * accel_ug_prtHz ./ (2*pi*freq_data).^2;
    
    lognoise = interp1(log10(freq_data),log10(disp_m_prtHz),log10(freq));
    noise_ASD = 10.^(lognoise);
    
    response_ZPK = tf(0);
    response_FR  = squeeze(freqresp(response_ZPK, 2*pi*freq));
    disp([' warning, the freq resp for the ' sensor_name ' is not yet defined']);
    
elseif strncmpi(sensor_name,'Wilcoxon_797L_spec',8)
    % https://buy.wilcoxon.com/amfile/file/download/file_id/85/product_id/429/
    % Meggitt/Wilcoxon model 797L
    % 100 mV/g = 0.010204 V / (m/s^2)
    % 2 Hz to 7 kHz
    
    freq_data = 10.^(1:3);
    accel_ug_prtHz = [8 5 5];
    disp_m_prtHz   = 9.8e-6 * accel_ug_prtHz ./ (2*pi*freq_data).^2;
    
    lognoise = interp1(log10(freq_data),log10(disp_m_prtHz),log10(freq));
    noise_ASD = 10.^(lognoise);
    
    response_ZPK = tf(0);
    response_FR  = squeeze(freqresp(response_ZPK, 2*pi*freq));
    disp([' warning, the freq resp for the ' sensor_name ' is not yet defined']);
    
elseif   strncmpi(sensor_name,'T360',4)
%     Freq = [3e-4 0.003  0.01 0.03 0.1  1  3.7  10 40];
%     ANdb = [-161  -187 -191.5 -194.5 -194.5  -189.5 -180  -169 -146];
    T360V1B = [0.0001    -146.9  % data from Nanometrics for current T360
               0.00032   -161.9
               0.001     -176.5
               0.0032    -186.8
               0.01      -192
               0.032     -194.5
               0.1       -194.5
               0.32      -193
               1         -189.4
               1.78      -186.2
               3.2       -182
               10        -168.5
               32        -150.5
               56        -141
               100       -129.5];
           
    Freq =   T360V1B(:,1);
    ANdb = T360V1B(:,2);
    
    Accnoisedb = interp1(log(Freq), ANdb,log(freq),'pchip');
    Accnoise_ASD = 10.^(Accnoisedb/20);
    noise_ASD = Accnoise_ASD./(2*pi*freq).^2;
    response_ZPK = tf(0);
    response_FR  = squeeze(freqresp(response_ZPK, 2*pi*freq));
    disp([' warning, the freq resp for the ' sensor_name ' is not yet defined']);
    
elseif   strncmpi(sensor_name,'T120H',5)
    T120H = [0.0001        -143.41
             0.00012115    -145.89
             0.00014678    -148.37
             0.00017783    -150.84
             0.00021544    -153.3
             0.00026102    -155.75
             0.00031623    -158.17
             0.00038312    -160.55
             0.00046416    -162.89
             0.00056234    -165.17
             0.00068129    -167.36
             0.0008254     -169.44
             0.001         -171.39
             0.0012115     -173.19
             0.0014678     -174.83
             0.0017783     -176.3
             0.0021544     -177.61
             0.0026102     -178.78
             0.0031623     -179.84
             0.0038312     -180.79
             0.0046416     -181.66
             0.0056234     -182.46
             0.0068129     -183.21
             0.008254      -183.91
             0.01          -184.56
             0.012115      -185.17
             0.014678      -185.73
             0.017783      -186.26
             0.021544      -186.73
             0.026102      -187.17
             0.031623      -187.55
             0.038312      -187.9
             0.046416      -188.2
             0.056234      -188.46
             0.068129      -188.67
             0.08254       -188.85
             0.1           -188.99
             0.12115       -189.09
             0.14678       -189.15
             0.17783       -189.18
             0.21544       -189.16
             0.26102       -189.11
             0.31623       -189.01
             0.38312       -188.87
             0.46416       -188.67
             0.56234       -188.41
             0.68129       -188.09
             0.8254        -187.68
             1             -187.19
             1.2115        -186.61
             1.4678	       -185.92
             1.7783        -185.12
             2.1544        -184.2
             2.6102        -183.15
             3.1623        -181.97
             3.8312        -180.65
             4.6416        -179.17
             5.6234        -177.51
             6.8129        -175.65
             8.254         -173.58
             10            -171.28
             12.115        -168.75
             14.678        -166.02
             17.783        -163.11
             21.544        -160.06
             26.102        -156.91
             31.623        -153.69
             38.312        -150.42
             46.416        -147.12
             56.234        -143.81
             68.129        -140.49
             82.54         -137.17
             100           -133.84
             121.15        -130.51
             146.78        -127.18
             177.83        -123.84
             215.44        -120.49
             261.02        -117.14
             316.23        -113.78
             383.12        -110.4
             464.16        -107
             562.34        -103.56
             681.29        -100.05
             825.4         -96.443
             1000          -92.657];
         
    Freq =   T120H(:,1);
    ANdb = T120H(:,2);
    
    Accnoisedb = interp1(log(Freq), ANdb,log(freq),'pchip');
    Accnoise_ASD = 10.^(Accnoisedb/20);
    noise_ASD = Accnoise_ASD./(2*pi*freq).^2;
    response_ZPK = tf(0);
    response_FR  = squeeze(freqresp(response_ZPK, 2*pi*freq));
    disp([' warning, the freq resp for the ' sensor_name ' is not yet defined']);
    
elseif   strncmpi(sensor_name,'T120C',5)
    T120C = [0.001   -145
             0.003   -154
             0.01   -161
             0.04   -168
             0.1   -171
             0.5   -174
             1   -175
             3   -173
             6   -169
             10   -166
             17.8   -160
             31.6    -153.5
             56.2    -144.5
             100   -134
             178    -124
             316    -114
             562    -103.5
             1000   -91];
         
    Freq =   T120C(:,1);
    ANdb = T120C(:,2);
    
    Accnoisedb = interp1(log(Freq), ANdb,log(freq),'pchip');
    Accnoise_ASD = 10.^(Accnoisedb/20);
    noise_ASD = Accnoise_ASD./(2*pi*freq).^2;
    response_ZPK = tf(0);
    response_FR  = squeeze(freqresp(response_ZPK, 2*pi*freq));
    disp([' warning, the freq resp for the ' sensor_name ' is not yet defined']);
    
elseif strncmpi(sensor_name,'BRS',3)
	% tilt noise fit based on BRS data - info added from SCModel.m 
    % points for f > 10 Hz are completely made up.
	freq_data  = [0.001  0.01   0.06   0.1    10      100      1000];  
	noise_data = [1e-8  1e-8  2e-10  1.5e-10  1.5e-10 1.5e-10  1.5e-10];

	lognoise =  interp1(log10(freq_data),log10(noise_data),log10(freq));
	noise_ASD   = 10.^lognoise;
    
    
    disp(['          warning, the output for ' sensor_name ' has units of rad/rtHz (not m/rtHz)']);
    
elseif strncmpi(sensor_name,'cBRS',4)
    % Pulled from SEI aLOG 1802, "cBRSNoiseModel.m"
    Q=1e2;
    f0=30e-3;

    lowangle=1e-11./(2*pi*freq).^(3);    
    bumpangle=2e-9./(1+freq./4e-1).^2;   
    trans=-(freq.^2)./(freq.^2-f0^2*(1+1i/Q));
    readoutNoise=abs(9e-11./trans);

    noise_ASD=sqrt(readoutNoise.^2+lowangle.^2+bumpangle.^2);
    
    disp(['          warning, the output for ' sensor_name ' has units of rad/rtHz (not m/rtHz)']);
    
elseif strncmpi(sensor_name,'CRS',3)
    % Pulled from https://dcc.ligo.org/LIGO-G2301261
    % page 8, "Instrument Noise" purple curve
    % (These are new results with HoQI readouts, better than Sep 2022
    % results shown in SEI aLOG 1967)
    % 1e-3, 1e2, and 1e3 data points completely made up
    % f > 2 Hz asymptote is a compromise between instrument noise and goal
    % assuming readout noise will be shot noise limited, but some fraction
    % worse than the goal for now.

    
    CRS = [ 1e-3        1.2e-06
            1e-2        1.2e-08
            1e-1        1.2e-10
            2e-1        4.0e-11
            5e-1        1.0e-11
            1e+0        6.0e-12
            2e+0        5.0e-12
            1e+1        5.0e-12
            1e+2        5.0e-12
            1e+3        5.0e-12];
        freq_data = CRS(:,1);
    disp_rad_prtHz = CRS(:,2);
             
    lognoise = interp1(log10(freq_data),log10(disp_rad_prtHz),log10(freq));
    noise_ASD = 10.^(lognoise);
    
    response_ZPK = tf(0);
    response_FR = 0;
    
    disp(['          warning, the output for ' sensor_name ' has units of rad/rtHz (not m/rtHz)']);
     
elseif  strncmpi(sensor_name,'PICO',4)
    % PICO PICOscale interfermoter from
    % http://www.smaract.com/wp-content/uploads/2017/01/PS-SS00007_ResolutionAndStability-1.pdf
    % see SEI aLOG 1366
    ss =  log10(5/4)/log10(0.009/0.01);
    Freq =  [0.001  0.009 0.01  0.04  0.1   0.3    1       10     100     1000];
    Disp  = [(0.009/0.001)^-ss*5e-9   5e-9 4e-9 2.3e-10 8.6e-11 5e-11 2.5e-11  1e-11  2.5e-12  1e-12 ];
    Dispnoise = interp1(log(Freq), log(Disp),log(freq),'pchip');
    noise_ASD = exp(Dispnoise);
    response_ZPK = tf(0);
    response_FR = 0;
    
elseif strncmpi(sensor_name,'SmarAct',7)
    % "Compact Michelson interferometers with subpicometer sensitivity"
    % Figure 2 (b), "Sensing head (FSS on)"
    % 1e-3 Hz and f > 30 Hz data points completely made up; low frequency
    % follows 1/f^3 slope, and high frequency assumes the 
    % sensor noise follows the laser amplitude noise.
    
    SmarAct = [1e-3      1.0e-08
               1e-2      1.0e-10
               2.5e-2    1.0e-11
               4e-2      4.0e-12
               1e-1      1.0e-12
               2e-1      4.0e-13
               3e-1      3.5e-13
               7e-1      3.0e-13
               1e+0      3.0e-13
               1e+1      3.0e-13
               1e+2      3.0e-13
               1e+3      3.0e-13];
           
    freq_data = SmarAct(:,1);
    disp_m_prtHz = SmarAct(:,2);
             
    lognoise = interp1(log10(freq_data),log10(disp_m_prtHz),log10(freq));
    noise_ASD = 10.^(lognoise);
    
    response_ZPK = tf(0);
    response_FR = 0;
    
    
elseif strncmpi(sensor_name,'HoQI',4)
    % 2018 performance from "A compact, large-range interferometer for
    % precision measurement and inertial sensing"
    % https://arxiv.org/pdf/1710.05943
    % Figure 3
    % 1e-3 and 1e3 Hz data point completely made up
    
    HoQI = [ 1e-3         1.0e-8 
             1e-2         1.0e-10
             3e-2         1.0e-11
             7e-2         1.0e-12
             1e-1         7.0e-13
             2e-1         3.0e-13
             1e+0         1.5e-13
             4e+0         1.0e-13
             1e+1         7.0e-14
             3e+1         4.0e-14
             7e+1         2.0e-14
             1e+2         2.0e-14
             1e+3         2.0e-14];
    
    freq_data = HoQI(:,1);
    disp_m_prtHz = HoQI(:,2);
             
    lognoise = interp1(log10(freq_data),log10(disp_m_prtHz),log10(freq));
    noise_ASD = 10.^(lognoise);
    
    response_ZPK = tf(0);
    response_FR = 0;

elseif strncmpi(sensor_name,'SPI_L',5)
    % Drawn from G2300846, page 17, Non-monolithic Interferometer, 50%
    % contrast
    % f > 1e1 Hz data points completely made up, assuming noise limit is 
    % shot-noise and thus flat, white
    
    SPI_L = [    1e-3                     1e-07
                 2e-3                   2.5e-08
                 1e-2                   1.5e-10
               1.2e-2                     1e-10
                 2e-2                     4e-11
                 3e-2                     3e-11
                 4e-2                     2e-11
                 1e-1                   1.3e-11
                 3e-1                     1e-11
                 1e+0                     1e-11
                 1e+1                     1e-11
                 1e+2                     1e-11
                 1e+3                     1e-11];
    
    freq_data = SPI_L(:,1);
    disp_m_prtHz = SPI_L(:,2);
             
    lognoise = interp1(log10(freq_data),log10(disp_m_prtHz),log10(freq));
    noise_ASD = 10.^(lognoise);
    
    response_ZPK = tf(0);
    response_FR = 0;

elseif strncmpi(sensor_name,'SPI_PY',5)
    % Drawn from G2300846, page 27, Delta RZ Interferometer, out-of-loop
    % Yaw measurement in rad/rtHz
    % Assumes lump in data between 0.12 and 0.40 Hz is some real signal,
    % not sensor noise
    % Assumes that P is the same as Y
    % 1e-3, and f > 1e1 Hz data points are completely made up
    SPI_PY = [   1e-3                     2.0e-05
                 1e-2                     3.0e-08
                 2e-2                     3.4e-09
                 3e-2                     1.7e-09
                 1e-1                     2.5e-10
               2.5e-1                     1.0e-10
                 4e-1                     7.0e-11
                 1e+0                     3.0e-11
                 2e+0                     1.7e-11
                 4e+0                     1.0e-11
                 1e+1                     1.0e-11
                 1e+2                     1.0e-11
                 1e+3                     1.0e-11];
             
    freq_data = SPI_PY(:,1);
    disp_rad_prtHz = SPI_PY(:,2);
             
    lognoise = interp1(log10(freq_data),log10(disp_rad_prtHz),log10(freq));
    noise_ASD = 10.^(lognoise);
    
    response_ZPK = tf(0);
    response_FR = 0;
    
    disp(['          warning, the output for ' sensor_name ' has units of rad/rtHz (not m/rtHz)']);
    
        
elseif strncmpi(sensor_name,'projections',4)
    % see T2100336
    temp = load('SEI_sensor_noise_projections'); % this file is in SEI's Common/MatlabTools
    noise_ASD = temp.sensor_noise_projections;
    % so, that's a little awkward, but what else to do?
    % eg  
    % cps_noise_1mm = SEI_sensor_noise('ADE_1mm', freq);
    % proj = SEI_sensor_noise('projections');
    % stg1_cps_x_noise = proj.bsc.stg1.cps.x * cps_noise_1mm;
    


else    
    disp(' error in SEI_sensor_noise  ');
    disp(['   sensor_name ' sensor_name ' is invalid'])
    disp(' defined sensors are: ')
    disp('          Displacement: ''ADE_1mm'', ''ADE_p25mm'', ''Kaman_IPSmeas'',''Fine_Hugh''')
    disp('          1Hz geophone: ''L4C'', ''GS13meas'', ''GS13calc'',');
    disp(' Low-freq Seismometers: ''T240meas'', ''T240spec'', ''Guralp_40T_Spec'',')
    disp('        Accelerometers: ''PCB_356M132_spec'', ''Wilcoxon_797L_spec'',')
    disp('      Rotation Sensors: ''BRS'', ''cBRS'', ''CRS'',')
    disp('       Interferometers: ''PICO'',''HoQI'',''SPI_L'',')
    disp('        Optical Levers: ''SPI_PY'',')
    disp('       Data Acquisiton: ''GS_16bit_ADC''.')
    disp(' still need to define ''STS2'' ')
    noise_ASD   = NaN * ones(size(freq));
    response_FR = NaN * ones(size(freq));
    
end


