function rt=GetTPRates(fn,phicol,psicol,outputfreq)
DHData = textread(fn);
size_DH = size(DHData);
Phis = DHData(:,phicol);
Psis = DHData(:,psicol);
s_s = zeros(size_DH(1),1);
j=1;
totS1=0;
totS2=0;
for i=1:size_DH(1),
    if Phis(i)>40 && Phis(i)<100
        break;
    else
        s_s(j)=TheState(Phis(i),Psis(i));
        j=j+1;
    end
end
state_space=s_s(1:(j-1));
size_DH(1)=max(size(state_space));
[X,pathsAB,pathsBA]=TransitionPaths(state_space);
%size of X gives us how many path I got.
%but I need to multiply using the frequency to 
%get the actual length in time.
size_x = size(X)
%pTP is the fraction of time I am on a transition
%along that entire trajectory.
%pTP = 0;
%for i=1:size_x(2),
%    pTP = pTP + X(i);
%end
pTP = sum(X)/size_DH(1)
%sum of the length of all paths
t_tp = mean(X)*outputfreq
t_tpsd = sqrt(var(X.*outputfreq))
max(X)*outputfreq
min(X)*outputfreq
X0=X;
X=X.*outputfreq;
rt=pTP/t_tp
%path_length_sum = 0;
%for i=1:size_x(2),
%    X(i) = X(i)*outputfreq;
%    path_length_sum = path_length_sum + X(i);
%end
%t_tp = path_length_sum/size_x(2)
%pTP = pTP/size_DH(1)
%####
% fid = fopen('states.dat','w');
% for i=1:size_DH(1)
%     fprintf(fid,'%2i\n',state_space(i));
% end
% fclose(fid);
% spab=size(pathsAB);
% fid1 = fopen('pathsAB.dat','w');
% for i=1:spab(1),
%    fprintf(fid1,'%i %i\n',pathsAB(i,1),pathsAB(i,2));
% end
% fclose(fid1);
% fid2 = fopen('pathsBA.dat','w');
% spba = size(pathsBA);
% for i=1:spba(1),
%     fprintf(fid2,'%i %i\n',pathsBA(i,2), pathsBA(i,1));
% end
% fclose(fid2);
%