Formula to calculate a t-distribution I am writing an application that will be dealing with <30 observations in a normal distribution. My understanding is that this point I would need to use t-distribution. The thing is, this is easy enough to look up in a table. However, I need this to be something to calculate it programmatically. I can't seem to find a formula to find the t-distribution score.
How can I build this table so to speak?  I am writing this in PHP.
I should also note that I am by no means a math/stat expert. So if my use of terms/vocabulary seems off just bear with me.
 A: My first impression is to use PHP build-in function link. I wonder why you cannot use the PHP stat extension?
If you really cannot use that extension, I suggestion you to generate a table and then use interpolation methods. For example, to get a cumulative distribution function table (CDF) to t-distribution, create a table like this t-distribution table.
A: Calculating and Plotting T-distribution 
You can use the following formulas to calculate t-distribution.
$ f(t)=\displaystyle\frac{1}{\sqrt{\nu}\,\mathrm{B}(\frac{1}{2},\frac{\nu}{2})}\left(1+\frac{t^2}{\nu}\right)^{-\frac{\nu+1}{2}} $
where:
$t$ = t value,
$\nu$ = degrees of freedom.
$\mathrm{B}(x,y)=\displaystyle\frac{(x-1)!(y-1)!}{(x+y-1)!}$
This is not a complex formula, you have to calculate factorial in the second formula which is usually written as $4!=4\times3\times2\times1$. This only valid if x and y are integers. 
Here I approach it in R:
t = 10
v = 3
x = 1/2
y = v/2 
B_x_y = (factorial((x-1)) * factorial((y-1)))/ factorial((x+y-1))
# when x and y are positive integers, otherwise we need more complex formula 
# we can just use beta function to derive this 
# B_x_y = beta(x,y)

# t-distribution 
f_t = (1/(sqrt(v)* B_x_y) *(1 + t^2/v)^(-((v+1)/2)))

# all of above can be done simply by using inbuilt R function dt 
# dt(10,3)

If you like to plot the distribution where df = 3. You can loop the process:
ft1 <- rep(NA, 100) 
tv <- seq(-6, 6, 0.1)

for ( i in 1:length(tv)){
  t = tv[i]
  v = 3
  x = 1/2
  y = v/2 
  B_x_y <- (factorial((x-1)) * factorial((y-1)))/ factorial((x+y-1))
  ft1[i] <- (1/(sqrt(v)* B_x_y) *(1 + t^2/v)^(-((v+1)/2)))
}         

plot(tv, ft1, type = "l", col = "blue", xlab = "t", ylab = "p", 
     main = paste ("t distribution with ", v, " degrees of freedom" )


All of the above process can be simply written as using dt function:
plot(tv, dt(tv, 3), col = "blue", type = "l", xlab = "t", ylab = "p") 

Creating Table
For creating table you need cumulative density function.  "The cumulative distribution function (CDF) describes the probability that a real-valued random variable X with a given probability distribution will be found to have a value less than or equal to x".
$ \displaystyle\int_{-\infty}^t f(u)\,du = \tfrac{1}{2} + t\frac{\Gamma \left( \tfrac{1}{2}(\nu+1) \right)} {\sqrt{\pi\nu}\,\Gamma \left(\tfrac{\nu}{2}\right)} \, {}_2F_1 \left( \tfrac{1}{2}, \tfrac{1}{2}(\nu+1); \tfrac{3}{2}; -\tfrac{t^2}{\nu} \right) $
${}_2F_1$ is a particular case of the hypergeometric function.
# pt(q, df)
# q is quantile
# Say we want 99% probability value (alpha = 0.01) for 3 df (looking table)
qt(0.99, 3)
4.540703

A: Easiest way is to use a software package to calculate this for you.
R is a standard platform for statistics and is free.
The main functions you would want to look up are pt, qt (You can get the manual of a function by using '?function' in R command line).
There are various wrappers to access R from other programming environments as well as other stats packages.
