clear all
close all

root = ['~/geospace_work/'];

station_list={'ott','iqa','frn'};


for station = station_list
  cd([root char(station)])
  disp(['Processing ' char(station)])
  load ([char(station) '_correlations_B_dB.mat']);

  X=[log10(maxB-minB); log10(maxBx-minBx); log10(maxBy-minBy)]';
  y= (maxdBdt);

  [i j]=find(isinf(X));
  X(i,:)=[]; y(i)=[];
  [i j]=find(isnan(X));
  X(i,:)=[];y(i)=[];
  
  f=find(log10(y)>4.9);
  y(f)=[];
  
  threshold_list=quantile(y/60,[0.6 0.7 0.8 0.9 0.95]);
  
  filename=[char(station) '_thresholds.mat']
  threshold_label={'0.6','0.7','0.8','0.9','0.95'};
  save(filename,'threshold_list','threshold_label');
  
end  
