yep = input('load data?');
if yep
  clear all
  disp('loading data');
  load('fp_max_real.dat');
  disp('Loaded fp_max_real');
  %  format [k,B,rho,mu,real,imag];
  load('rho.vdat');
  disp('Loaded rho');
  load('B.vdat');
  disp('Loaded B');
  load('k.vdat');
  disp('Loaded k');
end

disp('Using pre-loaded data');

%aalpha = 1; % f
%aalpha = -1; % d
%d = 1;

disp('opening file mu_max.dat');
fid = fopen('mu_max_log_2.dat','w');

disp('starting loop');
for i = 2:length(k)-1
   i                                   
   idx1 = find(abs(max_real(:,1) - k(i)) < sqrt(eps)); 
   Bvec = sort(max_real(idx1,2));
   Bvec = union(Bvec,Bvec);
   for j  = 1:length(Bvec)                                
   idx2 = find(abs(max_real(idx1,2) - Bvec(j)) < sqrt(eps));  
      for n = 1:length(rho)                              
         idx = find(abs(max_real(idx1(idx2),3) - rho(n)) < sqrt(eps));
         if isempty(idx); 
            disp(sprintf('Search empty for (%d,%d,%d)',i,j,n));
         end
         [val,midx] = max(max_real(idx1(idx2(idx)),5));
         fprintf(fid,'%24.16e',max_real(idx1(idx2(idx(midx))),:));
         fprintf(fid,'\n');                      
      end
   end
end
