function site=MakeLikelihood(site,deg,USS,datum,formatflag)


% where formatflag indicates:
% 1=years ad
% 2=years before datum
% 3=ka before datum
% 4=De scale

figure(1); clf;

%fit polynomials- Likelihoods stored on De scale
for i=1:length(site)
    [site(i).likelihood pairsUSS]=get_polyfit(site(i).pairs,deg,USS);
    
    % Plotting ----------------------------------------------
    subplot(length(site),1,i)
    %plot(site(i).pairs(:,1),site(i).pairs(:,2),'xk','markersize',4,'markeredgecolor',[0.5 0.5 0.5]);
    plot(pairsUSS(:,1),pairsUSS(:,2),'xk','markersize',4,'markeredgecolor',[0.5 0.5 0.5]);
    hold on
    plot(site(i).likelihood(:,1),site(i).likelihood(:,2),'-b')
    maxy=max(site(i).likelihood(:,2));
    axis([0 inf 0 2*maxy]);
    xlabel('De')
    ylabel(site(i).name)
    
end
clear priors


% Write the text file:
for i=1:length(site)
    
    t1=site(i).likelihood;
    t1(:,1)=t1(:,1)/site(i).dr;
    
    switch formatflag
        case 1
            t1(:,1)=datum-t1(:,1)*1000; % Gregorian year AD
        case 2
            t1(:,1)=t1(:,1)*1000;   % years before datum
        case 3
            t1(:,1)=t1(:,1);        %ka before datum
        case 4
            t1(:,1)=t1(:,1)*site(i).dr;
    end
            
            
    t1=sortrows(t1);
    step=t1(2,1)-t1(1,1);
    
    t2=[];
    t2(1,:)=[t1(1,1)-step 0];
    t2(2:length(t1(:,1))+1,:)=t1;
    t2(length(t1(:,1))+2,1:2)=[t1(length(t1(:,1)),1)+step 0];
    %priors(i).text=t2;
    
    
    textfile=[site(i).name,'.prior'];
    fid = fopen(textfile,'wt');
    for k=1:length(t2)
        fprintf(fid,'%4.8f %4.8f\n',t2(k,1),t2(k,2));
    end
        status = fclose(fid);
end

end



    


function [xycurve pairsUSS]=get_polyfit(pairs,deg,USS)
%xi is vector of points to interpolate fit
%deg is the degree of polynomial

pairs(:,1)=randn(length(pairs(:,1)),1).*(USS*pairs(:,1))+pairs(:,1);
pairs_sorted=sortrows(pairs);
x=pairs_sorted(:,1);
y=pairs_sorted(:,2);

xi=[min(x):(max(x)-min(x))/100:max(x)]';
p = polyfit(x,log(y),deg);
f = exp(polyval(p,x));

yi = interp1(x,f,xi,'linear');

ind1=find(yi>0);
ind2=not(isfinite(yi));
yi(ind2)=0;
xycurve=[xi(min(ind1):max(ind1)) yi(min(ind1):max(ind1))];
pairsUSS=pairs;
end