I'm trying to learn Quaternions for an implementation of L-Systems in 3D. I can't dumb it down using the pre-built Processing functions because I want to take the code and port it into Python as a module for Blender.

I have looked at the following resources and I think I understand axis and angle (self-taught in maths unfortunately) but I don't know how to effect a rotation, I'm not sure what my first and second Quaternions should be.  Each tutorial confuses me at the finish line by stating the solution as an equation I already don't understand.

http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/transforms/index.htm#sample

http://www.cprogramming.com/tutorial/3d/quaternions.html

The classes are from Ariel's ArcBall, with some of the math from the cprogramming link.

I would like to simply walk +100 pixels along the x axis and then rotate HALF_PI through the x axis (like being moved on the end of a stick). I would consider this a "Hello Quaternions!" I have posed this question with the classes below in the hope that someone might be able to show me how to apply them.
Code:

Vec3 here, there, axis;
float angle;
Quat rot;
void setup(){
  size(400, 400, P3D);
  here = new Vec3(0, 0, 0);
  there = new Vec3(100, 0, 0);
  axis = new Vec3(100, 0, 0);
  angle = HALF_PI;
  rot = new Quat();
}
void draw(){
  translate(200, 200);
  pushMatrix();
  translate(here.x, here.y, here.z);
  box(30, 30, 30);
  popMatrix();
  pushMatrix();
  translate(there.x, there.y, there.z);
  box(30, 30, 30);
  popMatrix();
}

static class Vec3{
  float x, y, z;
  Vec3(){
  }
  Vec3(float x, float y, float z){
    this.x = x;
    this.y = y;
    this.z = z;
  }
  void normalize(){
    float length = length();
    x /= length;
    y /= length;
    z /= length;
  }
  float length(){
    return (float) Math.sqrt(x * x + y * y + z * z);
  }
  static Vec3 cross(Vec3 v1, Vec3 v2){
    Vec3 res = new Vec3();
    res.x = v1.y * v2.z - v1.z * v2.y;
    res.y = v1.z * v2.x - v1.x * v2.z;
    res.z = v1.x * v2.y - v1.y * v2.x;
    return res;
  }
  static float dot(Vec3 v1, Vec3 v2){
    return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
  }
  static Vec3 mul(Vec3 v, float d){
    Vec3 res = new Vec3();
    res.x = v.x * d;
    res.y = v.y * d;
    res.z = v.z * d;
    return res;
  }
  void sub(Vec3 v1, Vec3 v2){
    x = v1.x - v2.x;
    y = v1.y - v2.y;
    z = v1.z - v2.z;
  }
}

static class Quat{
  float w, x, y, z;
  Quat(){
    reset();
  }
  Quat(float w, float x, float y, float z){
    this.w = w;
    this.x = x;
    this.y = y;
    this.z = z;
  }
  void reset(){
    w = 1.0f;
    x = 0.0f;
    y = 0.0f;
    z = 0.0f;
  }
  void set(float w, Vec3 v){
    this.w = w;
    x = v.x;
    y = v.y;
    z = v.z;
  }
  void set(Quat q){
    w = q.w;
    x = q.x;
    y = q.y;
    z = q.z;
  }
  static Quat mul(Quat q1, Quat q2){
    Quat res = new Quat();
    res.w = q1.w * q2.w - q1.x * q2.x - q1.y * q2.y - q1.z * q2.z;
    res.x = q1.w * q2.x + q1.x * q2.w + q1.y * q2.z - q1.z * q2.y;
    res.y = q1.w * q2.y + q1.y * q2.w + q1.z * q2.x - q1.x * q2.z;
    res.z = q1.w * q2.z + q1.z * q2.w + q1.x * q2.y - q1.y * q2.x;
    return res;
  }
  float[] getValue(){
    // transforming this quat into an angle and an axis vector...
    float[] res = new float[4];
    float sa = (float) Math.sqrt(1.0f - w * w);
    if (sa < EPSILON){
	sa = 1.0f;
    }
    res[0] = (float) Math.acos(w) * 2.0f;
    res[1] = x / sa;
    res[2] = y / sa;
    res[3] = z / sa;
    return res;
  }
}
Quat localRotation(float angle, Vec3 axis){
  return new Quat(
    cos(angle/2),
    axis.x * sin(angle/2),
    axis.y * sin(angle/2),
    axis.z * sin(angle/2));
}
float [] quatTo2Matrix(Quat q){
  float [] m = new float[16];
  m[0] = 1 - 2 * sq(q.y) - 2 * sq(q.z);
  m[1] = 2 * q.x * q.y - 2 * q.w * q.z;
  m[2] = 2 * q.x * q.z + 2 * q.w * q.y;
  m[3] = 0;
  m[4] = 2 * q.x * q.y + 2 * q.w * q.z;
  m[5] = 1 - 2 * sq(x) - 2 * sq(z);
  m[6] = 2 * q.y * q.z - 2 * q.w * q.z;
  m[7] = 0;
  m[8] = 2 * q.x * q.z - 2 * q.w * q.y;
  m[9] = 2 * q.y * q.z - 2 * q.w * q.z;
  m[10] = 1 - 2 * sq(x) - 2 * sq(y);
  m[11] = 0;
  m[12] = 0;
  m[13] = 0;
  m[14] = 0;
  m[15] = 1;
}
 

Can't say I can help you much, I'm on the same boat as you

So, I will cheerlead you on. This is something I wish to learn as well.

Ah, thanks. This is the trouble with developing code from two different websites. I'll have a play with that when I get home.

I started looking at the Matrix class that Processing uses. I noticed that transformation seems to work differently to the OpenGL site's method. It's a += of some multiplied matrix elements rather than add to the transform and then multiply by it.

The big trouble with my transformations is that each one multiplies on the last. Whereas in Processing, translate() does not apply any scale. I'm clearly doing it wrong.

Does anybody know what exactly I'm screwing up on?