%Find input parameter ranges for central case
%and Convert to natural units

clear all ; 

N=5000;
halfN=N/2;
density=15000; % 1/micron**2
h=0.05  ; % microns

%central case
TfCenter=1207 %microseconds
AmpCenter=873 %channels

serial=(1:1568)';

%=====The normal data=========================
fname='norm_avg_blip_results1.txt'
fullname = ( strcat('c:\My Documents\BigData\bartol\april2000runs\May22format\',fname) );

[str1,str2,tf,ba,R,kp,D,peak,tr20,tr80,tpeak,tf90,tf33] = ...  
   textread(fullname,...
   '%s%s%d%d%f%f%f%f%f%f%f%f%f',...
   1568);
%kp is in units of 1e7 1/Mole sec
%D is in units of 1e-6 cm**2/sec
%tf is in mSec

alpha = (1+R)./ (tf*.001);
beta = ba .* alpha;
km = beta ./(2*R);

%Tkm= 1e6 ./km ;
%Talphabeta = 1e6 ./(alpha+beta);

Tbind = 1.39 * h ./ (density * kp*1e7 * 1e15/6.02e23) *1e6;
Tdiff = 0.8/(4*pi) * N ./(density*D*1e-6*1e8) *1e6 ;

%symmetry = (tpeak - tr80)./(tf90 - tr80);
peakNorm = peak;
eff = peak ./ halfN ;
risetime = tr80 - tr20;
tr20norm=3*tr20./risetime;
falltime = tf33 - tf90;
plateautime = (tf90 - tr80)./sqrt(risetime.^2+falltime.^2);

%=====The DFP data=========================
fname='dfp_avg_blip_results.txt'
fullname = ( strcat('c:\My Documents\BigData\bartol\april2000runs\May22format\',fname) );

[str1,str2,tf,ba,R,kp,D,peak,tr20,tr80,tpeak,tf90,tf33] = ...  
   textread(fullname,...
   '%s%s%d%d%f%f%f%f%f%f%f%f%f',...
   1568);
peakDFP = peak;
effDFP = peak ./ halfN ;
risetimeDFP = tr80 - tr20;
falltimeDFP = tf33 - tf90;
tr20DFP=3*tr20./risetimeDFP;
plateautimeDFP = (tf90 - tr80)./sqrt(risetimeDFP.^2+falltimeDFP.^2);

%=====The BTX data=========================
fname='hbtx_avg_blip_results.txt'
fullname = ( strcat('c:\My Documents\BigData\bartol\april2000runs\May22format\',fname) );

[str1,str2,tf,ba,R,kp,D,peak,tr20,tr80,tpeak,tf90,tf33] = ...  
   textread(fullname,...
   '%s%s%d%d%f%f%f%f%f%f%f%f%f',...
   1568);

peakBTX=peak;
effBTX = peak ./ halfN ;
risetimeBTX = tr80 - tr20;
falltimeBTX = tf33 - tf90;
tr20BTX=3*tr20./risetimeBTX;
plateautimeBTX = (tf90 - tr80)./sqrt(risetimeBTX.^2+falltimeBTX.^2);

%=====Do conversions to real units===================

ampScale = peakNorm ./ AmpCenter;
timeScale = falltime ./ TfCenter;

Nreal = N ./ ampScale;

alpha = alpha .* timeScale;
beta  = beta .* timeScale;
Kp = kp .* timeScale;
Km = km .* timeScale;

Dreal = D .* timeScale./ampScale;

peakNorm = peakNorm./ampScale;
peakDFP = peakDFP./ampScale;
peakBTX = peakBTX./ampScale;

risetime = risetime ./ timeScale;
risetimeDFP = risetimeDFP ./ timeScale;
risetimeBTX = risetimeBTX ./ timeScale;

falltime = falltime ./ timeScale;
falltimeDFP = falltimeDFP ./ timeScale;
falltimeBTX = falltimeBTX ./ timeScale;

%=====Compute the acceptance ranges========

TargetRange = 0.2;  % +/- log units

for Target1=[2.6 ];%2.1 3.1]
   for Target2=[1.16 ];%1.70 0.70]
      for Target3=[0.69 ];%1.10 0.20]
         
         %=====The output data======================
         data = [ ba, R, kp, D,Tbind, Tdiff...
               beta, alpha, Kp, Km, Dreal, Nreal,...
               risetime,    .001*falltime,    peakNorm,plateautime,    tr20norm, eff, ...
               risetimeDFP, .001*falltimeDFP, peakDFP, plateautimeDFP, tr20DFP, effDFP, ...
               risetimeBTX, .001*falltimeBTX, peakBTX, plateautimeBTX, tr20BTX, effBTX, ...
               log(risetime ./ falltime), ...
               log(falltimeDFP ./ falltime), ...
               log(risetimeBTX ./ risetimeDFP) ...
               zeros(1568,1),... 
               serial,...
            ];
         
         
         %eliminate rows with eff<.05
         data = data(find(data(:,15)>=0.05*AmpCenter),:);  
         
         %get rows in Tr/Tf acceptance range
         data = data(find(-data(:,31)<(Target1+TargetRange)),:);
         data = data(find(-data(:,31)>(Target1-TargetRange)),:);
         
         %get rows in TfDFP/TfNorm acceptance range
         data = data(find(data(:,32)<(Target2+TargetRange)),:);
         data = data(find(data(:,32)>(Target2-TargetRange)),:);
         
         %get rows in TrBTX/TrDFP acceptance range
         data = data(find(data(:,33)<(Target3+TargetRange)),:);
         data = data(find(data(:,33)>(Target3-TargetRange)),:);
         
         %Compute sum(abs(log(target)))
         data(:,34) = sqrt( ...
            (data(:,31)+Target1).^2 ...
            + (data(:,32)-Target2).^2 ...
            +(data(:,33)-Target3).^2 );
         
         %sort by first 4 columns
         sortedData=sortrows(data,[34]);
         
         
         %normalize each parameter (column) to the range 0-1.0
         
         permute = [1 2 3 4];
         
         % 2 1 3 4         
         % 1 4 2 3
         
         for i=1:length(permute)
            pos(:,permute(i)) = .1*rand(size(data(:,i))) + ...
               (log(data(:,i))-log(min(data(:,i))))...
               /(log(max(data(:,i)))-log(min(data(:,i))));
         end
         pos = pos';
         
         npoints = length(data) ; 
         
         mapres = 10 ;
         map = jet(mapres) ; %define a color map jet
         
         dep=34; %index of the dependent variable to plot
         
         s = (mapres)*...
            (data(:,dep)-min(data(:,dep)))/((max(data(:,dep))-min(data(:,dep))))+1 ;
         %s = 1:npoints;
         highbrk = npoints;
         lowbrk = 1;
         %set(gcf,'position',[300 0 600 800] ) %make a bigger figure
         %orient tall   %fill the paper when we print
         
         cindex = fix(s(lowbrk));
         plot(pos(:,lowbrk), 'color', 'white' );
         
         %make the axes
         axisx = [1,4; 1,1; 2,2; 3,3; 4,4 ]';
         axisy = [0,0; 0,1; 0,1; 0,1; 0,1 ]';
         line(axisx, axisy, 'color', 'k','erasemode','none');
         drawnow;
         
         axisx=[1,2 ; 2,3; 3,4 ]';
         set(gca, 'drawmode','fast')
         
         for i = highbrk:-1:lowbrk
            if (s(i)<12)
               cindex = max(1,min(fix(s(i)),mapres)) ;
               axisy=[pos(1,i),pos(2,i); pos(2,i),pos(3,i);...
                     pos(3,i),pos(4,i)]';
               line(axisx,axisy,'erasemode','none','color',map(cindex,:),'linewidth',2);
               drawnow
            end
         end
         
      end
   end
end
%colorbar
title ('Color is distance to target -- blue=close')
set(gca,'xtick',[1 2 3 4])
set(gca,'xticklabel',['b/a|b/2Km|Kp|D'])
ylabel('normalized parameter range')