The model has been described fully elsewhere (International J. Mass and Heat Transfer, 43, 2059-2073 (2000). Authors are Harrison, Londry and Slavin.) The following is the documented SigmaPlot regression code.
 

alphae=.93                                                        'each "alpha" represents the corresponding
alphai=1                                                            'alpha in the theory.
alphas=2
alphar=1
'alphaer=.931                                                    'from P= 0 fit
adh=2.05e-5                                                     'equals alphac*delta / hr from P=0 fit
alphao=-1-ln((B/R)^2/2 + 2*lambda/R)
a=1.0                                                                'accommodation coefficient
sigma=gav/2.06                                                 'variance of the Gaussian distribution
pi=3.1416
D=2.18e-10                                                     'diameter of the He atom
k=1.38e-23                                                      'Boltzmann's constant
v=sqrt(8*k*T/(pi*4*1.66e-27))                        'average velocity of a gas molecule
R=0.0006                                                         'radius of a spheroid
eps=0.75                                                          'emissivity of the pebble
sig=5.67e-8                                                      'Stefan-Boltzmann constant

f1=0.0221                                                        'fraction of gaps in #1 of the 11 Gaussian bins
f2=0.0456
f3=0.0797
f4=0.1191
f5=0.1514
f6=0.164
f7=0.1514
f8=0.1191
f9=0.0797
f10=0.0456
f11=0.0221

Cv=1.5*k                                                         'heat capacity per gas molecule
n=P*1000/(k*T)                                               'number density of the gas molecules
lambda=1/(sqrt(2)*n*pi*D*D)                          'mean free path of a gas molecule
beta=(a/(3*(2-a)) - (1.5*25*pi/64))*n*v*Cv    'correction term for the gas conductivity
                                                                         'going from lambda << gav, to >> gav

Kg = (25*pi/64)*n*v*lambda*Cv                     'gas conductivity when lambda >> gav
Go=Kg*(1-exp(-R/lambda))*alphao*pi*R        'conductance for gas in the "outer gap"
Ks=3.33-3.088e-3*T+1.247e-6*T^2               'conductivity for bulk alumina (T=temp. in K)
Gs=Ks*alphas*0.5*pi*R                                   'conductance through the spheroid
Gc=Ks*adh/alphae                                            'conductance through the contacts
Gr=alphar*16*sig*R^2/(2/eps - 1)*T^3            'conductance by radiation between spheroids

'Distribute gaps into 11 Gaussian bins of width 2sigma/11, centred on gav.
 

g1=gav+2*sigma*(-5)/5                                     'width of the gap, for Gaussian bin #1
Kgprime1=Kg*(1-exp(-1.5*g1/lambda)) + beta*g1*exp(-1.5*g1/lambda)
                                                                          'conductivity of the gas in gap #1
L1=f1*Kgprime1/g1                                           'conductance from gap #1
 

g2=gav+2*sigma*(-4)/5
Kgprime2=Kg*(1-exp(-1.5*g2/lambda)) + beta*g2*exp(-1.5*g2/lambda)
L2=f2*Kgprime2/g2

g3=gav+2*sigma*(-3)/5
Kgprime3=Kg*(1-exp(-1.5*g3/lambda)) + beta*g3*exp(-1.5*g3/lambda)
L3=f3*Kgprime3/g3

g4=gav+2*sigma*(-2)/5
Kgprime4=Kg*(1-exp(-1.5*g4/lambda)) + beta*g4*exp(-1.5*g4/lambda)
L4=f4*Kgprime4/g4

g5=gav+2*sigma*(-1)/5
Kgprime5=Kg*(1-exp(-1.5*g5/lambda)) + beta*g5*exp(-1.5*g5/lambda)
L5=f5*Kgprime5/g5

g6=gav+2*sigma*0/5
Kgprime6=Kg*(1-exp(-1.5*g6/lambda)) + beta*g6*exp(-1.5*g6/lambda)
L6=f6*Kgprime6/g6

g7=gav+2*sigma*1/5
Kgprime7=Kg*(1-exp(-1.5*g7/lambda)) + beta*g7*exp(-1.5*g7/lambda)
L7=f7*Kgprime7/g7

g8=gav+2*sigma*2/5
Kgprime8=Kg*(1-exp(-1.5*g8/lambda)) + beta*g8*exp(-1.5*g8/lambda)
L8=f8*Kgprime8/g8

g9=gav+2*sigma*3/5
Kgprime9=Kg*(1-exp(-1.5*g9/lambda)) + beta*g9*exp(-1.5*g9/lambda)
L9=f9*Kgprime9/g9

g10=gav+2*sigma*4/5
Kgprime10=Kg*(1-exp(-1.5*g10/lambda)) + beta*g10*exp(-1.5*g10/lambda)
L10=f10*Kgprime10/g10

g11=gav+2*sigma*5/5
Kgprime11=Kg*(1-exp(-1.5*g11/lambda)) + beta*g11*exp(-1.5*g11/lambda)
L11=f11*Kgprime11/g11
 

Gi=(L1+L2+L3+L4+L5+L6+L7+L8+L9+L10+L11)*alphai*pi*B^2      'total conductance
                                                                                                                'from all gaps
Gc=0                                                                         'Gc not included; too small
Gpar=Gi +Gc
f=(0.5*alphae/R)*(Gr+Go + Gs*Gpar/(Gs+Gpar))     'Overall function to fit to experimental data.
fit f to K                                                                      'SigmaPlot command to fit f to the expt'l
                                                                                   'thermal conductivity data labelled K.

It is also necessary to define the columns which represent K, P and T, and provide starting parameters for gav and B.