function out=flyLIF2(params,connect,inputs,inputPSCs,TRPCurrent)


numNeurons=size(connect,1);
numInputs=length(inputs);


% TRPCurrent=zeros(78,1);
% TRPCurrent(50:59)=pA*10^-9;

TRPCurrent=TRPCurrent*10^-9;

spikeThr=zeros(numNeurons,1)-45e-3;

if numInputs ~=size(inputPSCs,2)
    error('input current count mismatch');
end

if isfield(params,'time');
    time=params.time; else time=0.5; end;                   % 0.5 (or 0.083s) seconds
if isfield(params,'dt');
    dt=params.dt; else dt=1e-4; end;                        % 1e-4 s
if isfield(params, 'N');
    N=params.N; else N=time/dt; end;
if isfield(params,'V0');
    V0=params.V0; else V0=-52e-3; end;                      % -52e-3 (or -40e-3, -60e-3) V
if isfield(params,'spikeMin');
    spikeMin=params.spikeMin; else spikeMin=-72e-3; end;    %-72e-3 (or -90e-3) V
% if isfield(params,'spikeThr');
%     spikeThr=params.spikeThr; else spikeThr=-45e-3; end;    %-45e-3 (or -30e-3) V
if isfield(params,'Cmem');
    Cmem=params.Cmem; else Cmem=2e-8; end;                  %2e-8F
if isfield(params,'Rmem');
    Rmem=params.Rmem; else Rmem=1e6; end;                   %1e6 Ohm
if isfield(params,'apDuration');
    apDuration=params.apDuration; else apDuration=2; end;   %2 (or 0.5, 8) ms
if isfield(params,'pscRiseT');
    pscRiseT=params.pscRiseT; else pscRiseT=0.002; end;     %2 ms
if isfield(params,'pscFallT');
    pscFallT=params.pscFallT; else pscFallT=0.005; end;                 %5 ms
if isfield(params,'pscMag');
    pscMag=params.pscMag; else pscMag=5e-9; end;            %5e-9 (or 2.5e-9, 20e-9)A

if isfield(params,'pscWeights');
    pscWeights=params.pscWeights; else pscWeights=1; end;            %5e-9 (or 2.5e-9, 20e-9)A
if length(pscWeights)==1
    inputWeights=zeros(numInputs,1)+pscWeights;
end


spikeMax=20e-3;     %purely cosmetic parameter

apRiseTime=round((apDuration*10)/2);
apRise=normpdf(-1:1/apRiseTime:0);
apRise=(apRise-min(apRise))/(max(apRise)-min(apRise));
apRise=apRise*(spikeMax-spikeThr(1))+spikeThr(1);
apFallTime=round((apDuration*9)/2);
apFall=sin(pi()/2:pi()/apFallTime:3*pi()/2);
apFall=(apFall-min(apFall))/(max(apFall)-min(apFall));
apFall=apFall*(spikeMax-spikeMin)+spikeMin-.0001;
AP=[apRise apFall]';
AP(1)=[];
AP(end+1)=AP(end)+0.0001;
AP(end+1)=AP(end)+0.0001;
apLength=length(AP);

pscRise=sin(linspace(-pi()/2,pi()/2,pscRiseT/dt));
pscRise=(pscRise-min(pscRise))/(max(pscRise)-min(pscRise));
pscFall=2.^(-(0:(7*pscFallT/dt))*dt/pscFallT);
pscFall=(pscFall-min(pscFall))/(max(pscFall)-min(pscFall));
PSC=[pscRise pscFall]';
PSC=PSC*23.6305/sum(PSC); %normalize PSC magnitude to counteract changes in duration or time-constants
PSC=PSC*pscMag/max(PSC);

tBuff=length(PSC);

%prepare input current vectors
Iin=zeros(N+tBuff,numNeurons);
for i=1:numInputs
    if ~ iscell(inputPSCs)
        if min(inputPSCs(:))==0
        pscFallTimes=find(inputPSCs(:,i)==1);
        else
            pscFallTimes=inputPSCs(:,i);
        end
    else
        pscFallTimes=inputPSCs{i};
    end
    for j=1:length(pscFallTimes)
        PSCrange=pscFallTimes(j):(pscFallTimes(j)+tBuff-1);
        Iin(PSCrange,inputs(i))=Iin(PSCrange,inputs(i))+PSC*inputWeights(i);
    end
end

V=zeros(N+tBuff,numNeurons)+V0;
V(1,:)=V0;
I=zeros(N+tBuff,numNeurons);
APMask=zeros(N+tBuff,numNeurons);

for t=2:N
    for i=1:numNeurons
        if APMask(t,i)==0
            if V(t-1,i)>spikeThr(i)
                V(t-1,i)=spikeThr(i);
                V(t:t+apLength-1,i)=AP;
                APMask(t:t+apLength-1,i)=1;
                I(t:t+tBuff-1,i)=I(t:t+tBuff-1,i)+PSC;
            else
                currentTemp=Iin(t-1,i) + sum(connect(:,i).*I(t-1,:)') + TRPCurrent(i);
                dV=(1/Cmem)*( (V0-V(t-1,i))/Rmem + currentTemp );
                V(t,i)=V(t-1,i)+dV*dt;
            end
        end
    end
end

out.Iin=Iin;
out.V=V;
out.I=I;
out.inputPSCs=inputPSCs;
out.APMask=APMask;
out.connect=connect;
out.P=params;
out.P.time=time;
out.P.dt=dt;
out.P.V0=V0;
out.P.spikeMin=spikeMin;
out.P.spikeThr=spikeThr;
out.P.Cmem=Cmem;
out.P.Rmem=Rmem;
out.P.apDuration=apDuration;
out.P.pscRiseT=pscRiseT;
out.P.pscFallT=pscFallT;