% Nernst
% with numberMolculesIn & movement based on charge Gradient
% With Ion Potential
% Determine the grid size and generate the diffusion.
  N = input('Enter Time:');
  nRedMolecules=20;  %Number of Potassium Particles
  permRed=0.375;  %permeability of Potassium ions
  nBlueMolecules=20; %sodium
  permBlue=0.2; %
  nGreenMolecules=20;
  permGreen=0.1; %permeability of Chloride ions
  membrane=0; %Membrane at x=0
%K
  posRed='In';
  [a b]=DiffusionMemPos(N,membrane,permRed,posRed);
  xRed=[a];
  xRedLength=length(xRed);
  yRed=[b];
  redIn=18; %number K inside cell initially
  for red=2:nRedMolecules
      if red>redIn-0.5;
          posRed='Out';
      end
      [c d] = DiffusionMemPos(N,membrane,permRed,posRed);
      xRed=[xRed;c];
      yRed=[yRed;d];
  end
  %Track of molecules in and out
  [positiveInR positiveOutR] =numberMoleculesIn(N,membrane,xRed,permRed,nRedMolecules);
%Na
  posBlue='In';
  [e f]=DiffusionMemPos(N,membrane,permBlue,posBlue); %Initial array for Na particles
  xBlue=[e]; %Initial X matrix for Na Particles
  xBlueLength=length(xBlue);
  yBlue=[f]; %Initial y Martix for Na particles
  blueIn=2; %Number of Na inside cell initially
  for blue=2:nBlueMolecules
      if blue>blueIn-0.5
          posBlue='Out';
      end
      [g h] = DiffusionMemPos(N,membrane,permBlue,posBlue);
      xBlue=[xBlue;g]; %Concatinate each x row
      yBlue=[yBlue;h]; %Concatinate each y row
  end
  %Track of molecules in and out
  [positiveInB positiveOutB] =numberMoleculesIn(N,membrane,xBlue,permBlue,nBlueMolecules);
  positiveIn=positiveInR + positiveInB;
  positiveOut=positiveOutR + positiveOutB;
  %Cl
  posGreen='In';
  [i j]=DiffusionMemPos(N,membrane,permGreen,posGreen); %Initial array for Cl particles
  xGreen=[i]; %Initial X matrix for Cl Particles
  xGreenLength=length(xGreen);
  yGreen=[j]; %Initial y Martix for Cl particles
  greenIn=4; %Number of Cl in initially
  for green=2:nGreenMolecules
      if green>greenIn-0.5
          posGreen='Out';
      end
      [k l] = DiffusionMemPos(N,membrane,permGreen,posGreen);
      xGreen=[xGreen;k]; %Concatinate each x row
      yGreen=[yGreen;l]; %Concatinate each y row
  end
  [negativeIn negativeOut] =numberMoleculesIn(N,membrane,xGreen,permGreen,nGreenMolecules);
  %Change Based On Charge Gradient & Membrane potential
  %Potassium
  [VK]=IonEquilibriumPotential(positiveInR, positiveOutR);
  [xRed1 chargeInR chargeOutR]=moveParticlesChargePot(N,membrane,xRed,permRed,nRedMolecules,positiveIn, positiveOut, negativeIn, negativeOut,VK);
  positiveIn=positiveIn + chargeInR;
  positiveInR=positiveInR + chargeInR;
  positiveOut=positiveOut + chargeOutR;
  positiveOutR=positiveOutR + chargeOutR;
  xRed=xRed1;
  %Sodium
  [VNa]=IonEquilibriumPotential(positiveInB, positiveOutB);
  [xBlue1 chargeInB chargeOutB]=moveParticlesChargePot(N,membrane,xBlue,permBlue,nBlueMolecules,positiveIn, positiveOut, negativeIn, negativeOut,VNa);
  positiveIn=positiveIn + chargeInB;
  positiveInB=positiveInB + chargeInB;
  positiveOut=positiveOut + chargeOutB;
  positiveOutB=positiveOutB + chargeOutB;
  xBlue=xBlue1;
  %Chloride
  [VCl]=IonEquilibriumPotential(negativeOut, negativeIn);
  [xGreen1 chargeInG chargeOutG]=moveParticlesChargePot(N,membrane,xGreen,permGreen,nGreenMolecules,positiveIn, positiveOut, negativeIn, negativeOut,VCl);
  negativeIn=negativeIn + chargeInG;
  negativeOut=negativeOut + chargeOutG;
  xGreen=xGreen1;
  %Ending Membrane Potentials
  [VK]=IonEquilibriumPotential(positiveInR, positiveOutR);
  [VNa]=IonEquilibriumPotential(positiveInB, positiveOutB);
  [VCl]=IonEquilibriumPotential(negativeOut, negativeIn);
  % Create the figure window
  close all
  figure(1)
  set(gcf,'position',[150 50 600 600])
  hold on
  
% Boundary
M = 20;
% Set the axes
  axis([-20 20 -20 20])
  axis equal square manual
% Label in and out
  yText=M-1;
  xTextMin=-M+.5;
  xTextMax=M-7.5;
  text(xTextMin,yText,'Inside Cell')
  text(xTextMax,yText,'Outside Cell')
% Set membrane
  xMem=membrane;
  yMin=-20;
  yMax=20;
% Animates the ions
  mov=avifile('Final2.avi');
  for col=2:xRedLength %By column
      plot([xMem xMem],[yMin yMax],'-c')
      for row=1:nRedMolecules %By Row
          plot(xRed(row,col-1),yRed(row,col-1),'.w','Markersize',20)
          plot(xRed(row,col),yRed(row,col),'.r','Markersize',20)
      end
      for row=1:nBlueMolecules %By Row
          plot(xBlue(row,col-1),yBlue(row,col-1),'.w','Markersize',20)
          plot(xBlue(row,col),yBlue(row,col),'.b','Markersize',20)
      end
      for row=1:nGreenMolecules %By Row      
          plot(xGreen(row,col-1),yGreen(row,col-1),'.w','Markersize',20)
          plot(xGreen(row,col),yGreen(row,col),'.g','Markersize',20)
      end
      pause(.1)
      F=getframe(gcf);
      mov=addframe(mov,F);
  end
  mov=close(mov);
  %Plot Membrane Potentials
  figure(2);clf;
  hold on
  xVMin=0;
  xVMax=2;
  Vm=58*log((permRed*(positiveOutR)+permBlue*(positiveOutB)+permGreen*(negativeIn))/(permRed*(positiveInR)+permBlue*(positiveInB)+permGreen*(negativeOut)));
  plot([xVMin xVMax],[Vm Vm],'m')
  text(1,Vm,'Vm')
  plot([xVMin xVMax],[VK VK],'r')
  text(1,VK,'VK')
  plot([xVMin xVMax],[VNa VNa],'b')
  text(1,VNa,'VNa')
  plot([xVMin xVMax],[VCl VCl],'g')
  text(1,VCl,'VCl')
  ylabel('Voltage(mV)')