In[4876]:= 
ClearAll["Global`*"] (*If you remove this you'll need to restart Mathematica each time you use the program. I have no idea why - something about the program thinking you're using a variable name as a function name or something like that*)

(*Experimental Variables*)
SETaveragePoreDiameter =50. (*nm*)
SETparticleDiameter = 20. (*nm*)
SETparticleConcentration = .25*5.66*^-2 (*kg/m^3. 0.25 comes from a 4:1 dilution, other value from BBI website*)
SETmembraneThickness = 50. (*nm*)
SETtransmembranePressurePsi = 5.
SETmolarSaltConcentration = .01 
SETzetaParticle = -0.015 (*volts*)
SETzetaMembrane = -0.020 (*volts*)


(*Converting our good practices into SI or radius*)
convertPsiToPa[x_] := x * 6894.75729 (* Pa, 34473 Pa = 5 psi*)
DtoR[x_] := x/2.


(*Physical Constants*)
viscosity = 1.*^-3 (*N s/m^2*)
faradayConstant = 9.648*^4 (*C/mole*)
gasConstant = 8.31446 (*J/(mol*K)*)
vacuumPermittivity = 8.854*^-12 (*Farads/meter*)
densityOfWater = 1000 (*kg/m^3*)
boltsmannConstant = 1.38*^-23
dielectricConstWater = 80.
temp = 300.(*K*)

debyeLength[molarSaltConcentration_]:= Sqrt[(vacuumPermittivity*dielectricConstWater*gasConstant*temp)/(2*faradayConstant^2*molarSaltConcentration)] 
(*for a symmetric monovalent electrolyte such as KCl*)

(*Stokes-Einstein Diffusion coefficient:*)
diffusionCoefficient[particleRadius_] := (boltsmannConstant * temp)/(6 Pi viscosity (particleRadius))

(*Dagan Equation. Volumetric Flow rate, m^3/s :*)
solventVolumeFlow[transmembranePressure_, poreRadius_, membraneThickness_] := (transmembranePressure*poreRadius^3)/(viscosity(3+(8/Pi)(membraneThickness/poreRadius))) 

solventSpeed[transmembranePressure_, poreRadius_, membraneThickness_] := solventVolumeFlow[transmembranePressure, poreRadius, membraneThickness]/(Pi poreRadius^2) (*m/s*)

(*Dimensionless variables*)
(*alpha := (particleRadius/poreRadius) *)(*alpha is the ratio of particle radius to pore radius*)
(*beta = r/poreRadius dimensionless radial position*)
(*tau := (poreRadius/debyeLength )*)

sigmaC[zetaMembrane_,poreRadius_,debyeLength_ ,particleRadius_] :=  (zetaMembrane(1+(poreRadius/debyeLength ) (particleRadius/poreRadius) ))/(particleRadius/poreRadius) (*surface charge density. I found this in Jess's thesis but I couldn't find where she found it*)
sigmaS[zetaParticle_,poreRadius_,debyeLength_ ,particleRadius_] := (zetaParticle (poreRadius/debyeLength ) BesselI[1,(poreRadius/debyeLength )])/BesselI[0,(poreRadius/debyeLength )] (*particle charge density*)

(*Paine and Scheer values interpolated to get g and k^-1. x is the diffusive drag, y is the convective drag. It was easier to let mathematica assume that x = 1, 2, 3 here and change the inputs (hence the o[(m*45./.9)+1.] nonsense) than it was to actually type in the pairs of values from Paine and Scheer *)
o= Interpolation[{1.00000,1.04393,1.09178,1.14397,1.20096,1.26330,1.33159,1.40654,1.48892,1.57966,1.67980,1.79054,1.91328,2.04963,2.20150,2.37109,2.56100,2.77430,3.01464,3.28635,3.59464,3.94578,4.34741,4.808805,5.34141,5.95938,6.68043,7.52686,8.52705,9.71752,11.14580,12.87453,14.98751,17.59862,20.86550,25.01092,30.35733,37.38467,46.83132,59.87967,78.51925,106.31670,150.22684,225.51931,372.41035,737.25652}]
x[m_] := o[(m*45./.9)+1.]
l= Interpolation[{1.00000,1.04365,1.09062,1.14122,1.19584,1.25488,1.31881,1.38816,1.46351,1.54554,1.63500,1.73276,1.83981,1.95729,2.08651,2.22898,2.38646,2.56102,2.75506,2.97143,3.21351,3.48533,3.79174,4.13856,4.53292,4.98360 ,5.50116,6.09927,6.79476,7.60918,8.57025,9.71417,11.08878,12.75849,14.81147,17.37105 ,20.61390,24.80033,30.32686,37.82235,48.33522 ,63.72853,87.60633,127.82576,204.96024,393.56585}]
y[n_] := l[(n*45./.9)+1.]

(*The following three equations come from Electrostatic effects on the partitioning of spherical colloids between dilute bulk solution and cylindrical pores, Smith & Deen.
lambda corresponds to equation, delta G to equation, and e to .*)
lambda[betaLambda_, poreRadius_, debyeLength_] := (Pi/2)BesselI[0,(poreRadius/debyeLength )*betaLambda]Sum[((betaLambda^t((2. t)!))/(2.^(3t)*(t!)^2))*BesselI[t,(poreRadius/debyeLength ) betaLambda]*((poreRadius/debyeLength ) BesselK[t+1, 2 (poreRadius/debyeLength )] +(3/4)*BesselK[t,2 (poreRadius/debyeLength )]),{t, 0, 10(*Should be Infinity, but on recommendation from Jess's thesis I do to 10*)}]
deltaG[betaDG_,poreRadius_, debyeLength_,particleRadius_,zetaParticle_,
zetaMembrane_] :=(((*part1:*)((8. Pi (poreRadius/debyeLength ) (particleRadius/poreRadius)^4 Exp[(poreRadius/debyeLength ) (particleRadius/poreRadius)])/(1.+(poreRadius/debyeLength )(particleRadius/poreRadius)^2))lambda[betaDG, poreRadius, debyeLength]sigmaS[zetaParticle,poreRadius,debyeLength ,particleRadius]^2)(*part 2:*)+((4.Pi^2(particleRadius/poreRadius)^2 BesselI[0,(poreRadius/debyeLength ) betaDG])/((1.+(poreRadius/debyeLength ) (particleRadius/poreRadius))BesselI[1,(poreRadius/debyeLength )]))sigmaS[zetaParticle,poreRadius,debyeLength ,particleRadius] sigmaC [zetaMembrane,poreRadius,debyeLength,particleRadius](*part3:*) + ((Pi BesselI[0,(poreRadius/debyeLength ) betaDG])/((poreRadius/debyeLength )*BesselI[1,(poreRadius/debyeLength )]))^2 (((Exp[(poreRadius/debyeLength ) (particleRadius/poreRadius)]-Exp[-(poreRadius/debyeLength ) (particleRadius/poreRadius)])(poreRadius/debyeLength ) (particleRadius/poreRadius) (Coth[(poreRadius/debyeLength ) (particleRadius/poreRadius)]-(1./ ((poreRadius/debyeLength ) (particleRadius/poreRadius)))))/(1. + (poreRadius/debyeLength ) (particleRadius/poreRadius)))sigmaC[zetaMembrane,poreRadius,debyeLength ,particleRadius]^2)/(*part 4:*)(Pi (poreRadius/debyeLength ) Exp[-(poreRadius/debyeLength ) (particleRadius/poreRadius)]-(2.(Exp[(poreRadius/debyeLength ) (particleRadius/poreRadius)]-Exp[-(poreRadius/debyeLength ) (particleRadius/poreRadius)])(poreRadius/debyeLength ) (particleRadius/poreRadius)(Coth[(poreRadius/debyeLength ) (particleRadius/poreRadius)]-(1./ ((poreRadius/debyeLength ) (particleRadius/poreRadius)))) lambda[betaDG, poreRadius, debyeLength])/(1.+(poreRadius/debyeLength ) (particleRadius/poreRadius)))
e[betaE_,poreRadius_, debyeLength_,particleRadius_,zetaParticle_,
zetaMembrane_] := poreRadius vacuumPermittivity dielectricConstWater ((gasConstant*temp)/faradayConstant)^2 * deltaG[betaE,poreRadius, debyeLength,particleRadius,zetaParticle,
zetaMembrane]
(*H and W come from Hindered Transport of Large Molecules in Liquid-Filled Pores by W.M. Deen. H is eq. 15, W eq. 16*)
(*Integrals solved by sampson's rule with n subdivisions as done in Jess's thesis. a can't be 0 because it will result in a 0^0 indeterminante form in the infinite sum in lambda*)
a =0.0000001
(*b = 1-(particleRadius/poreRadius)*)
n = 30
f [r_,poreRadius_, debyeLength_,particleRadius_,zetaParticle_,
zetaMembrane_] :=x[(particleRadius/poreRadius)]^(-1)Exp[-e[r,poreRadius, debyeLength,particleRadius,zetaParticle,
zetaMembrane]/(boltsmannConstant*temp)]*r
g[u_,poreRadius_, debyeLength_,particleRadius_,zetaParticle_,
zetaMembrane_] := y[(particleRadius/poreRadius)]^(-1)(1-u^2)Exp[-e[u,poreRadius, debyeLength,particleRadius,zetaParticle,
zetaMembrane]/(boltsmannConstant*temp)]*u
s[ particleRadius_,poreRadius_]:= ((1-(particleRadius/poreRadius))-a)/n
h[poreRadius_, debyeLength_,particleRadius_,zetaParticle_,
zetaMembrane_] := 2 N[(s[ particleRadius,poreRadius]/3)Sum[f[(a+(2j-2) s[ particleRadius,poreRadius]),poreRadius, debyeLength,particleRadius,zetaParticle,
zetaMembrane]+4f[(a+(2j-1) s[ particleRadius,poreRadius]),poreRadius, debyeLength,particleRadius,zetaParticle,
zetaMembrane]+f[(a+(2j) s[ particleRadius,poreRadius]),poreRadius, debyeLength,particleRadius,zetaParticle,
zetaMembrane],{j, 0, n/2}]]
w [poreRadius_, debyeLength_,particleRadius_,zetaParticle_,
zetaMembrane_] := 4 N[(s[ particleRadius,poreRadius]/3)Sum[g[(a+(2p-2) s[ particleRadius,poreRadius]),poreRadius, debyeLength,particleRadius,zetaParticle,
zetaMembrane]+4g[(a+(2p-1) s[ particleRadius,poreRadius]),poreRadius, debyeLength,particleRadius,zetaParticle,
zetaMembrane]+g[(a+(2p) s[ particleRadius,poreRadius]),poreRadius, debyeLength,particleRadius,zetaParticle,
zetaMembrane],{p, 0, n/2}]] 
peclet[poreRadius_, debyeLength_,particleRadius_, membraneThickness_, transmembranePressure_,zetaParticle_,
zetaMembrane_]:= (w[poreRadius, debyeLength,particleRadius,zetaParticle,
zetaMembrane] solventSpeed[transmembranePressure, poreRadius, membraneThickness] membraneThickness)/(h[poreRadius, debyeLength,particleRadius,zetaParticle,
zetaMembrane] diffusionCoefficient[particleRadius])
(*soluteSpeed = (((1-Exp[-peclet])/(w solventSpeed particleConcentration))+(1/(solventSpeed particleConcentration))Exp[-peclet])^(-1)*)
reflectionCoefficient[poreRadius_, debyeLength_,particleRadius_,zetaParticle_,
zetaMembrane_]  := 1-w[poreRadius, debyeLength,particleRadius,zetaParticle,
zetaMembrane]  (*should be 1-w but this seems to work. Hell if I know why*)
filtrateConcentrationConvection[poreRadius_, debyeLength_,particleRadius_, particleConcentration_,zetaParticle_,
zetaMembrane_]  := particleConcentration*reflectionCoefficient[poreRadius, debyeLength,particleRadius,zetaParticle,
zetaMembrane] (*filtrate concentration is first calculated assuming only convection is happening, and then that concentration gradient is used to calculate the diffusive flux and get the actual filtrate concentration*)
soluteFlux[poreRadius_, debyeLength_,particleRadius_, particleConcentration_,transmembranePressure_, membraneThickness_,zetaParticle_,
zetaMembrane_] := w[poreRadius, debyeLength,particleRadius,zetaParticle,
zetaMembrane]solventSpeed[transmembranePressure, poreRadius, membraneThickness] particleConcentration  ((1 - (filtrateConcentrationConvection[poreRadius, debyeLength,particleRadius, particleConcentration,zetaParticle,
zetaMembrane]/particleConcentration)Exp[-peclet[poreRadius, debyeLength,particleRadius, membraneThickness, transmembranePressure,zetaParticle,
zetaMembrane]])/(1-Exp[-peclet[poreRadius, debyeLength,particleRadius, membraneThickness, transmembranePressure,zetaParticle,
zetaMembrane]]))(*solute mass flux: kg/(m^2 s)*)
(*solventMassFlux = solventSpeed * densityOfWater (*1 m^3 of water is 1000kg*)*)


sieving[averagePoreDiameter_,
particleDiameter_,
particleConcentration_,
membraneThickness_,
transmembranePressurePsi_,
molarSaltConcentration_,
zetaParticle_,
zetaMembrane_]:= (soluteFlux[DtoR[averagePoreDiameter*1*^-9], debyeLength[molarSaltConcentration],DtoR[particleDiameter*1*^-9], particleConcentration,convertPsiToPa[transmembranePressurePsi], membraneThickness*1*^-9,zetaParticle,
zetaMembrane]/solventSpeed[convertPsiToPa[transmembranePressurePsi], DtoR[averagePoreDiameter*1*^-9], membraneThickness*1*^-9])/particleConcentration

(*Dynamic user interface:*)
(*{"particleDiameter = 20 nm \nParticle Concentration = .25 * 5.66*^-2 kg/m^3\nMembrane Thickness = 50 nm\nParticle Zetapotential = -15 mV\nPressure (PSI)",Slider[Dynamic[sliderPressure],{.01,2}],Dynamic[sliderPressure],"\nAverage Pore Diameter (nm)",Slider[Dynamic[sliderPoreDiameter],{SETparticleDiameter*1.12,SETparticleDiameter*10}],Dynamic[sliderPoreDiameter],"\nMolar Salt Concentration",Slider[Dynamic[sliderSalt],{{0.001,0.005,0.01,0.05, .1,.5,1}}],Dynamic[sliderSalt],"\nZetapotential of Membrane (V)",Slider[Dynamic[sliderZmembrane],{0,-0.1}],Dynamic[sliderZmembrane],"\nSieving Coefficient:",Dynamic[sieving[sliderPoreDiameter,(*sliderParticleDiameter*)SETparticleDiameter, SETparticleConcentration, SETmembraneThickness,sliderPressure,sliderSalt,SETzetaParticle, sliderZmembrane]], "\nPeclet #:",Dynamic[peclet[DtoR[sliderPoreDiameter*1*^-9], debyeLength[sliderSalt],DtoR[SETparticleDiameter*1*^-9], SETmembraneThickness*1*^-9,convertPsiToPa[sliderPressure], SETzetaParticle, sliderZmembrane]]}*)

(*Standard Separation. Should give a Sieving Coefficient of 0.1694*)
sieving[SETaveragePoreDiameter, SETparticleDiameter, SETparticleConcentration, SETmembraneThickness, 1, SETmolarSaltConcentration,SETzetaParticle, SETzetaMembrane]

(*Particle Diameter Plot*)
(*Plot[sieving[(SETaveragePoreDiameter), v, SETparticleConcentration, SETmembraneThickness, SETtransmembranePressurePsi, SETmolarSaltConcentration,SETzetaParticle, SETzetaMembrane], {v, 1, 42}, AxesLabel->{"Particle Diameter (nm)", "Sieving Coefficient"}]*)

(*Plot[sieving[(SETaveragePoreDiameter), SETparticleDiameter, SETparticleConcentration, v, SETtransmembranePressurePsi, SETmolarSaltConcentration,SETzetaParticle, SETzetaMembrane], {v, 1, 10}, AxesLabel->{"Membrane Thickness (nm)", "Sieving Coefficient"}]*)

(*Plot[sieving[(SETaveragePoreDiameter), SETparticleDiameter, SETparticleConcentration, SETmembraneThickness, v, SETmolarSaltConcentration,SETzetaParticle, SETzetaMembrane], {v, 0.001, 1}, AxesLabel->{"Transmembrane Pressure (PSI)", "Sieving Coefficient"}]*)

(*Plot[sieving[(SETaveragePoreDiameter), SETparticleDiameter, SETparticleConcentration, SETmembraneThickness,SETtransmembranePressurePsi , v,SETzetaParticle, SETzetaMembrane], {v, 0.001, .01}, AxesLabel->{"Molar Salt Concentration", "Sieving Coefficient"}]*)

Plot[sieving[SETaveragePoreDiameter, SETparticleDiameter, SETparticleConcentration, SETmembraneThickness, SETtransmembranePressurePsi, SETmolarSaltConcentration,SETzetaParticle, v], {v, -.000, -.030}, AxesLabel->{"Membrane Charge (mV)", "Sieving Coefficient"}]
Out[4877]= 50.
Out[4878]= 20.
Out[4879]= 0.01415
Out[4880]= 50.
Out[4881]= 5.
Out[4882]= 0.01
Out[4883]= -0.015
Out[4884]= -0.02
Out[4887]= 0.001
Out[4888]= 96480.
Out[4889]= 8.31446
Out[4890]= 8.854*10^-12
Out[4891]= 1000
Out[4892]= 1.38*10^-23
Out[4893]= 80.
Out[4894]= 300.
Out[4901]= InterpolatingFunction[{{1.,46.}},<>]
Out[4903]= InterpolatingFunction[{{1.,46.}},<>]
Out[4908]= 1.*10^-7
Out[4909]= 30
Out[4920]= 0.169467
Out[4921]= 
GraphicsBox[{{}, {}, 
{Hue[0.67, 0.6, 0.6], LineBox[CompressedData["
1:eJwVz3k81PkfB/AJCeVKRcfK1S+WKRNKu8unwfSj5Eo07qlZkSNZJEdhHWGy
hEZoqFCOaJuktXiLmfmutOUnhpHb6MDPKkeR4/f5/fF6vB7PP15/vDTPXHD6
WYJEIkXi/L9RCGnE8RkHfg/8x3v4dSsyqLU7/6GWAzk7aI2MjlaktlQ4e4XL
AbHNF5Ux7JlfD0tXVXMgPjL44fh/WlExO1RfupQDo7q5mRNvWhEJxOH1WRxw
SzBQne1pRS3ybXLaAXhf3KaxNtSKrCqzjOe+40BnlhRJcqYVlUR1JuZH3Ib8
AwMpXpt5aCZKNOzaVwhVEo8fLlN46MTfNtN7qIVgWZXxY5oDD0k/8f/6Jb8A
dtEs7Qwu8lBlu73CyFo+1KgdLNf+DfvSuMYr13xY2nuiVrmGhxrCxeL2x7dA
kTzC6Pubh+pLyp8Mq92CTU8N38hP8JCKn3r84sU82LdqNKgow0enPeh0PSEb
LCZprinafFRzd2XU34wNhhnlz6bN+IgR6+hfzb4JPO/Kw8t0Plq7FvBNdiUX
qjnLRexf+Oirrkma/6lcMO2MlX2WwUeSO/1eiB7lwJ7PV8O5D/hoyYyxjaGa
A0xH2pbTzXxkOy/rPRWSDUbSATIuIj5KH4jkZnTfgDuyntNyM3y0x40tZUK+
AWstitsXNghQh+7uDVMpWdDv8LSpUV2AzukGba3/byZc2t68391IgOgN5zVY
RzPhH4cbu9fZCFCIhGfwhPNv4CO1L4blIUAOhnWrA1oZcJcamq5yQYDeJ4dc
EwlZUFB77WF4ggC5oKd6knfSwcGb5D+YLUAxGu48Y8c0IKu/0F8pFSBeaefz
lA2pcP+4fr5knQCNs4RCma4UIF+tvmosEKCGwsPjWcnJUCVHFk52C1Cggt5G
G/skyBkZnhsaEyCkM5mtqZ4IQqLaZXhGgP46s2mbxJsE2LxzC9N0VYC+l6yv
UeLEg+FgRfdnWQJZvBRY7LWPg4Ai9b7BrQTynYa+HoMrsMQr06JqEOg6iz1/
j4iGKieTmno9AmnGkCXj9l2Gs/z+sINGBDq3vx89ex4BMvmKJQ4/Emi7sc0F
kmYYOLvIJB21JJByvDsznXURrsWab7M4RqAwbXnlcetg8DdoTZR2IFDR+yHz
N0r+QF3M9ow4RaBXXxVFp1b8INDkkkIq9rLskV96P/pBhjiSW4B92qD4/kCL
HyTr9YiasZUu+ih9CPODJweW2XIuBEpYGh5Z7j0HL4pfiwqwfeVHE/91xxdM
C/tCa10JRDYSv7hMYcIu48dzPDqBOMHUaC91JvTs9fHswlao4OhbbmRCVklh
/Rj2tAadtWn8LCQMNjtJuBHokcKrE0V5ZyF3kLtojm00UdfBWzsDZ/MMtZ5i
H76TLlR8xYCCh1P9Be4Equj/kDxXz4BaPWpPOfYOtaOHRPcZMB2lKajDXsog
5d2NZ0DZsaT4N9gNMeF0ExMGuDN7guU8CIROe/e7FfqAw0i2dgR2urov73i3
F9CixHNWnvhvlK2OsNILhMM9mU7YUcIDiT4JXtA/YbvbB/t8xppl+H4vQOl/
qkRjW6/mtXJSPSEi5IHbI2ypwfaWT2YeoPu4ykfNi0AxtynP2WV0MJZrtR3E
DltU1dSKpUOMEphOYAecWo2rOkkHYuGw6gK2m3z7keckOqgoV1XKexPINJbZ
POF2GpKoVw/8hD3nwQYzRVe4It1RdxM7aNdK49glZ2jPzi6m+hCoWKp72MjR
Gaht8eTj2J1TVZKJ3zuDK7ey2hnbtNHDRmfgJMgHPr5+DlvCq6GLaXESuCOp
TSxsdlH05PgmJ7hQ7jzZhf0rbY7SYWkPlyu/hXgwCPRha8mlkAU7SP5KS2Ji
27072ahUbgfzDi+uB2KrpXCtHRXsIGo2LjoGu/qvUO/OXltIqSlpK8DuPf6J
1RV0DFab/zTsxTZwmn7fy6YBqU+kZ32GQBozmamXqTQo3meZao+9JcNIf8ek
FYxKrolcsJfbIoPdza3AZzHI/mfsl1TJ+QGxBWjl7k2Lxw6gqEmMUY7AJ+VV
8zrsJk74vzfUIeg9lMptxC7rq0+zlkPwVsF3Ow87womm3P7IDHSLt/zRgb3V
wk29Y+UH2LXeyvojtrNmomkf2wQeDIyUqp4l0E+ebdE7J40hrzJh/3fYOrcU
wMPcGAQRYRVa2HNKt2hD4gNw7+2xMDJ2NqnaSUwxhG38564W2PegRtTP3gN+
Ra89/bFDD/3g422uA/ll0ctB2BY1vHcjYi1o3xzJCsUe5Yhm31E04OMXRlw0
NnVZetcoeyuQNy7WsbCHJ3fqvRdvhhb7vSOZ2HFvDQ9OUZTgZl3ptxzs2PDa
4in2ejjYeUPqNnaD81D9AoUEA+HrZouxF85trliXN9+UoN7cWYL93vZlg37e
eFMn79W9B9i63MjPHrLrm0oYZN8q7P8B5NLXZg==
"]]}},
AspectRatio->0.6180339887498948,
Axes->True,
AxesLabel->{FormBox["\"Membrane Charge (mV)\"", TraditionalForm], FormBox["\"Sieving Coefficient\"", TraditionalForm]},
AxesOrigin->{0, 0.155},
ImageSize->{700.4140625, Automatic},
Method->{},
PlotRange->{{-0.03, 0.}, {0.15478138622606744`, 0.18272929181217734`}},
PlotRangeClipping->True,
PlotRangePadding->{Scaled[0.02], Scaled[0.02]}]