function [newsp,immunemn,newmn,newvar] = killingfn(x,T,mn,lim,dp,specmn,specvar,dpmx,klngspred,klngparam,maxind)

% The closer you are to the original mean, mn1, the more likely you will
% die...

%note that here we are including two different modes of resistance
%development(the two tails)

%these correspond to SNPs and copy number changes
%note however, that we could model this in more than one dimension
%(make the normal values not a vector but a matrix)
% and this would allow for more dimensions of disease resistance
% development

%SHOULD I INCLUDE TIME... I GUESS SO!

%variables: counter (p)checkdist (p)klngspred (p)klngparam
%(p)knockdownprobdeath (p)nrmcst dumdum (p)immunemn (p)grwrate (p)clonepop (p)carrycap (p)totalpop (p)sumcomb (p)newmn (p)stnd (p)newvar (p)newsp

%let's save knockdownprobdeath as a vector... so we can view it later!

knockdowndeath = zeros(length(x),1);

checkdist =  zeros(length(x),1);

ctr = 1;

for i = 1:length(x);
    
    checkdist(i)  = abs(x(i) - mn);
    % the knockdown is quite large, and then there is an immune system
    % suppression (that is, the change the survivors of the drug will survive)
    % is a negative exponential function of time
     
    %now apply the drugpressure:    
    
    knockdownprobdeath(i) = normpdf(checkdist(i),mn,klngspred) * (dp/dpmx);
    
    
    %so this is a spread out version of the original model, will allow
    %the setting of a value just above the value of the pdf
    %at the mean to work for entire distribution; note we use the mean
    %for the drug
    
    %Let's kill two standard deviations from mean... 
    %Chance of survival then is about .04 at mean
    %That is, the probabilty that they will die near the mean
    %is the probabilty that they are two standard devs from the mean
    
    %If they are two standard devs from the mean, the prob they will
    %survive is 1/sqrt(2*pi)
    
    nrmcst = (normpdf(0,mn,klngspred))/erf(klngparam/sqrt(2));
    
    %erf(klngparam/sqrt(2)) = the probability that the value is klngparam
    %standard deviations from the mean (mn); by dividng, we are killing
    %a certain percentage (e.g. 95.5%) at the mean
    
    
    
    %*********************
    %the closer you are to the drug mean(fixed), the more you are are punished
    %*********************
   
    %mindist = normpdf(0,specmn,specvar)
    %drugpressweight = dp*(nrmcst - mindist)
    
    if knockdownprobdeath(i) >= random('unif',0,nrmcst,1);
        
        x(i) = NaN;
        
    end;

%remove the dead
    
    if ~isnan(x(i));
        
        dumdum(ctr) = x(i); %dumdum is a row vector
        ctr = ctr + 1; 
    end
    
end;

%hist(x)
%title('This is a hist of the resulting species.')

%figure

%hist(dumdum)
%title('This is a hist with the dead removed.')

%test = waitforbuttonpress;

%if test == 0;
%reset counter

ctr = 1;

%Now immune system activates, kills off an additional number:
  
%Now we need to reduce the poplation by regrowing in competition:
%endemic condition so same as before...

%need to check if we should change the mean (mutants are more fit); 
%then drift mean

%random walk the mean

immunemn = mn;

if rand>=.99; 
    
    immunemn = immunemn+1; 

end

    if rand>=.99; 
        
        immunemn = immunemn-1; 
    end


%note that the immune suppression is based on the fact that those closer to
%the mean benefit...

% but what about the chance that a mutant is more fit?

%grow and cut

grwrate = zeros(length(dumdum),1);

for i = 1:length(dumdum);
    %*********************
    %the farther away you are from the mean, the more you are are punished
    %*********************
    
grwrate(i) = abs(dumdum(i) - immunemn);

end




clonepop=ones(1,length(dumdum)); %a row vector
carrycap = lim;
totalpop = 0;
fastermnmult = 2;

for t = 1:T
    
    for i = 1:length(dumdum);

        %grow them all
        
        clonepop(i) = round((1 + fastermnmult * grwrate(i)) * clonepop(i)); 
        totalpop = totalpop + clonepop(i);
        
    end

    if totalpop>=lim;
        
        
        for i = 1:length(dumdum);
            clonepop(i) = round(clonepop(i)/totalpop *lim);
        end
        
        totalpop=lim;
        
    end;
    
    for i = 1:length(dumdum);
        
    grwrate(i) = abs(grwrate(i)*(1 - totalpop/carrycap));
    
    end    
  
end


%clonepop;
%test = waitforbuttonpress;
%if test == 0;
%hist(clonepop);
%end;



%now we weight:

sumcomb = 0;

for i = 1:length(dumdum);
      
sumcomb = sumcomb + clonepop(i)*dumdum(i); 

end


newmn = sumcomb/length(dumdum);

stnd = std(clonepop.*dumdum);
newvar = stnd^2;


%resample
newsp = random('norm',newmn,newvar,1,maxind);


%variables:
%newsp bigdum dumdum counter i expcst nrmcst knockdownprobdeath mn 

clear dumdum counter i

%end