OK here is the shader source code:

GLcharARB* lorentz_vert_source =

“void main(void)”

“{”

" gl_TexCoord[0] = gl_MultiTexCoord0;"

" gl_TexCoord[1] = gl_MultiTexCoord1;"

" gl_TexCoord[2] = gl_MultiTexCoord2;"

" gl_TexCoord[3] = gl_MultiTexCoord3;"

" gl_Position = ftransform();"

“}”;

GLcharARB* lorentz_frag_source =

“uniform float time;”

“uniform float time_step;”

“uniform float Bo;”

“uniform float Re;”

“uniform float xscale;”

“uniform sampler2D velocity_mass;”

“uniform sampler2D position_charge;”

“uniform sampler2D vel_mass;”

“uniform sampler2D pos_charg;”

“float ERRTOL =0.08;”

“float ERRT =0.05;”

“uniform int particle_tex_height;”

“uniform int particle_tex_width;”

"

“//1

“vec3 ElectricField(vec3 position, int x){”

" const float one_over_four_pi_eo = 0.28;”

" int i, j;"

" vec2 particle_index;"

" vec4 particle_position;"

" vec3 direction;"

" float magnitude;"

" vec3 electric_sum = vec3(0.0);"

" for(i=0; i<int(particle_tex_height); i++) {"

" for(j=0; j<int(particle_tex_width); j++) {"

" particle_index.x = (float(j)+0.5) / float(particle_tex_width);"

" particle_index.y = (float(i)+0.5) / float(particle_tex_height);"

“if(x==0){”

" particle_position = texture2D(position_charge, particle_index);"

“} if(x==1){”

“particle_position = texture2D(pos_charg, particle_index);”

“}”

" direction = position - particle_position.xyz;"

" magnitude = length(direction);"

" if (magnitude > 1e-9) {"

" electric_sum += one_over_four_pi_eo * particle_position.w * direction / pow(magnitude, 2.0);"

" }"

" }"

" }"

" return electric_sum;"

“}”

"

" //2

“float EllipticFirstKind(in float x, in float y, in float z)”

“{”

“float alamb, ave, delx, dely, delz, sqrtx, sqrty,sqrtz,xt,yt,zt, ea, eb;”

“xt=x;”

“yt=y;”

“zt=z;”

“do {”

“sqrtx = sqrt(xt);”

“sqrty = sqrt(yt);”

“sqrtz = sqrt(zt);”

“alamb = sqrtx*(sqrty+sqrtz)+sqrty*sqrtz;"*

"xt = 0.25(xt+alamb);”

“yt = 0.25*(yt+alamb);”

“zt = 0.25*(zt+alamb);”

“ave = (1.0/3.0)*(xt+yt+zt);"*

“delx = (ave-xt)/ave;”

“dely = (ave-yt)/ave;”

“delz = (ave-zt)/ave;”

“} while (max(max(abs(delx),abs(dely)),abs(delz)) > ERRTOL);”

"ea=delxdely-delz*delz;"*

"eb=delxdely*delz;”

“return(1.0+((1.0/24.0)*ea-(0.1)-(3.0/44.0)*eb)**ea+(1.0/14.0)**eb)/sqrt(ave);"*

“}”

"

" //3

“float LEFK (in float ak)”

“{”

“float phi = 3.1415/2.0;”

“float s = sin(phi);”

"return sEllipticFirstKind(cos(phi)*cos(phi),(1.0-s*ak)(1.0+sak),1.0);”

“}”

"

" //4

“float EllipticSecondKind(in float x, in float y, in float z)”

“{”

```
"float alamb, ave, delx, dely, delz, sqrtx, sqrty,sqrtz,xt,yt,zt,fac,sum, ea, eb,ec,ed,ee;"
"xt=x;"
"yt=y;"
"zt=z;"
"sum=0.0;"
"fac=1.0;"
"do {"
"sqrtx = sqrt(xt);"
"sqrty = sqrt(yt);"
"sqrtz = sqrt(zt);"
"alamb = sqrtx*(sqrty+sqrtz)+sqrty*sqrtz;"
"sum += fac/(sqrtz*(zt+alamb));"
"fac=0.25*fac;"
"xt = 0.25*(xt+alamb);"
"yt = 0.25*(yt+alamb);"
"zt = 0.25*(zt+alamb);"
"ave = 0.2*(xt+yt+3.0*zt);"
"delx = (ave-xt)/ave;"
"dely = (ave-yt)/ave;"
"delz = (ave-zt)/ave;"
"} while (max(max(abs(delx),abs(dely)),abs(delz)) > ERRT);"
"ea = delx*dely;"
"eb = delz*delz;"
"ec = ea-eb;"
"ed = ea-6.0*eb;"
"ee = ed+ec+ec;"
"return 3.0*sum+fac*(1.0+ed*(-(3.0/14.0)+(0.25*(9.0/22.0))*ed-(1.5*(3.0/26.0))*delz*ee)+delz*((1.0/6.0)*ee+delz*(-(9.0/22.0)*ec+delz*(3.0/26.0)*ea)))/(ave*sqrt(ave));"
```

“}”

"

" //5

////Preforms Legendre Transformation to K space of Elliptic Integral Second Kind////////////////////////////

“float LESK (in float ak)”

“{”

“float cc,q,phi;”

“phi = 3.1415/2.0;”

“float s = sin(phi);”

“cc = cos(phi)*cos(phi);"*

"q = (1.0-sak)*(1.0+s*ak);”

“return s*(EllipticFirstKind(cc,q,1.0)-((s*ak)*(s*ak))*EllipticSecondKind(cc,q,1.0)/3.0);”

“}”

"

" //6

//CALCULATES MODULUS - SEE SITE ON OFF AXIS FIELD STRENGTH FOR EQUATION

“float Kv( float r, float a, float z)”

“{”

“return sqrt((4.0*r*a)/((z*z)+(r*r)+(a*a)+(2.0*r*a)));"*

“}”

"

" //7

“vec3 MagneticFieldRing(float x, float y, float z, float a, float I, float magz_pos)”

“{”

“float brr;”

“float bzz;”

“float posx = x;”

“float posy = y;”

“float posz = z-magz_pos;”

"float r = sqrt(((posxposx)+(posy*posy)));”

“float theta = atan(posy/posx);”

```
"float ak = Kv(r,a,posz);"
"brr = (I*posz/(2.*3.1415*r*sqrt((((a*a)+(r*r)+(2.*r*a))+posz*posz))))*(-LEFK(ak)+((a*a+r*r+posz*posz)/(((a*a)+(r*r)-(2.*a*r))+posz*posz))*LESK(ak));"
"bzz = ((I)/(2.*3.1415*sqrt((a*a)+(r*r)+(2.*r*a)+(posz*posz))))*(LEFK(ak)+((a*a-r*r-posz*posz)/(((a*a)+(r*r)-(2.*a*r))+posz*posz))*LESK(ak));"
"float b_x = brr*cos(theta);"
"float b_y = brr*sin(theta);"
"float b_z = bzz;"
"return vec3(b_x, b_y, b_z);"
```

“}”

"

" //8

“uniform int number_of_magnets;”

“uniform float magnet_current[6];”

“uniform float magnet_zpos[6];”

“uniform float magnet_radius[6];”

“vec3 MagneticField(float time, vec3 position)”

“{”

//DOD - this function is not correct

“float posx = position.x;”

“float posy = position.y;”

“float posz = position.z;”

“int j;”

//“float b_mag;”

```
"vec3 b_temp = vec3(0.0);"
"vec3 b = vec3(0.0);"
"for (j=0; j < 3; j++)"
"{"
"b = MagneticFieldRing(posx, posy, posz,magnet_radius[j],magnet_current[j],magnet_zpos[j]);"
"b_temp = b_temp + b;"
```

“}”

“return b_temp;”

“}”

"

" //10

“vec3 Lorentz(in float charge, in float mass, in float time, in vec3 position, in vec3 velocity, in vec3 electric_field)”

“{”

" float chrg_div_mass = charge/mass;"

" vec3 fields = electric_field + cross(velocity, MagneticField(time, position));"

" return chrg_div_mass * fields;"

“}”

"

" // 10

“vec3 Euler(vec3 y, vec3 slope, float delta)”

“{”

" return y + slope * delta;"

“}”

"

" //11

"

"

“void RungeKutta(float charge, float mass, float time, float time_step, in vec3 position, in vec3 velocity, in vec3 electric_field, out vec4 computed_vel, out vec4 computed_pos)”

“{”

```
" float half_step = 0.5 * time_step;"
" float half_time = time + half_step;"
" float full_time = time + time_step;"
" vec3 k1, k2, k3, k4;"
" vec3 ks1, ks2, ks3, ks4;"
" vec3 v1, v2, v3, v4;"
" vec3 p1, p2, p3, p4;"
" v1 = velocity;"
" p1 = position;"
" k1 = Lorentz(charge, mass, time, p1, v1, electric_field);"
" ks1 = v1*xscale;"
" v2 = Euler(velocity, k1, half_step);"
" p2 = Euler(position, ks1, half_step);"
" k2 = Lorentz(charge, mass, half_time, p2, v2,electric_field);"
" ks2 = v2*xscale;"
" v3 = Euler(velocity, k2, half_step);"
" p3 = Euler(position, ks2, half_step);"
" k3 = Lorentz(charge, mass, half_time, p3, v3,electric_field);"
" ks3 = v3*xscale;"
" v4 = Euler(velocity, k3, time_step);"
" p4 = Euler(position, ks3, time_step);"
" k4 = Lorentz(charge, mass, full_time, p4, v4,electric_field);"
" ks4 = v4*xscale;"
" computed_vel.xyz = velocity + (time_step / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4); "
" computed_vel.w = mass;"
" computed_pos.xyz = position + (time_step / 6.0) * (ks1 + 2.0 * ks2 + 2.0 * ks3 + ks4); "
" computed_pos.w = charge;"
```

“}”

"

" //12

“void main(void)”

“{”

" int x;"

" vec4 computed_velocity[2];"

" vec4 computed_position[2];"

" vec4 space[2];"

" vec4 momentum[2];"

" float mass[2];"

" float charge[2];"

" momentum[0] = texture2D(velocity_mass, gl_TexCoord[0].xy);"

" space[0] = texture2D(position_charge, gl_TexCoord[1].xy);"

" momentum[1] = texture2D(vel_mass, gl_TexCoord[2].xy);"

" space[1] = texture2D(pos_charg, gl_TexCoord[3].xy);"

" mass[0] = momentum[0].a;"

" mass[1] = momentum[1].a;"

" charge[0] = space[0].a;"

" charge[1] = space[1].a;"

" for(x=0;x<2;x++){"

“vec3 electric_field = ElectricField(space[x].xyz, x);”

“RungeKutta(charge[x], mass[x], time, time_step, space[x].xyz, momentum[x].xyz, electric_field, computed_velocity[x], computed_position[x]);”

“}”

//LOOP REDUCES TEMP REGISTER USE FROM 57 TO 41, ACCOUNTING FOR ALL 3 MAGNETS

" gl_FragData[0]= computed_velocity[0];"

" gl_FragData[1] = computed_position[0];"

" gl_FragData[2] = computed_velocity[1];"

" gl_FragData[3] = computed_velocity[1];"

“}”;

Anybody have any suggestions? Thank you