// spherical_deriv.java // given U(r,t,p) radius, theta in xy plane, phi from vertical axis // x = r sin p cos t // y = r sin p sin t // z = r cos p // // r = sqrt(x^2 + y^2 + z^2) r > 0 // t = atan(y/x) t = -Pi .. Pi t = 0 along positive x axis // p = acos(z/r) p = 0 ..Pi p = 0 along positive z axis // p = Pi along negative z axis // if p=0 or p=Pi then x=0 and y=0, thus force t=0 (tt) // // ^ ^ ^ // compute gradient: dU/dr r 1/r*sin(p) dU/dt t 1/r dU/dp p // // dx = sin p cos t dr - r sin p sin t dt + r cos p cos t dp // dy = sin p sin t dr + r sin p cos t dt + r cos p sin t dp // dz = cos p dr - r sin p dp // // dr = sin p cos t dx + sin p sin t dy + cos p dz // dt = -(cos t / r sin p) dx + (cos t/ r sin p) dy + 0 dz // dp = (cos p sin t /r) dx + (cos p sin t /r) dy - sin p /r dz // // test on simple outward normal ? = 1.0 // U(r,t,p) = r + r*t + r*sin(t)*p public class spherical_deriv { int nr = 16; int nt = 16; int np = 16; double rv[] = new double[nr]; // radius values at i double tv[] = new double[nt]; // theta values at j double pv[] = new double[np]; // phi values at k double U[][][] = new double[nr][nt][np]; // function values double Pi = Math.PI; public spherical_deriv() { double r, t, tt, p, x, y, z, rb, tb, pb, xb, yb, zb, err, err1; double r0 = 0.25; double dr = 2.0/(double)nr; double t0 = -Pi; // -Pi to Pi double dt = 2.0*Pi/(double)nt; double p0 = 0.0; double dp = Pi/(double)np; System.out.println("spherical_deriv.java running"); System.out.println("check functions"); x = 1.5; y = 2.2; z = 3.1; System.out.println("x="+x+", y="+y+", z="+z); r = fr(x,y,z); t = ft(x,y,z); p = fp(x,y,z); System.out.println("r="+r+", t="+t+", p="+p); xb = fx(r,t,p); yb = fy(r,t,p); zb = fz(r,t,p); err = Math.abs(x-xb)+Math.abs(y-yb)+Math.abs(z-zb); System.out.println("xb="+xb+", yb="+yb+", zb="+zb+", err="+err); rb = fr(xb,yb,zb); tb = ft(xb,yb,zb); pb = fp(xb,yb,zb); err = Math.abs(r-rb)+Math.abs(t-tb)+Math.abs(p-pb); System.out.println("rb="+rb+", tb="+tb+", pb="+pb+", err="+err); System.out.println(" "); x = -1.5; y = 2.2; z = -3.1; System.out.println("x="+x+", y="+y+", z="+z); r = fr(x,y,z); t = ft(x,y,z); p = fp(x,y,z); System.out.println("r="+r+", t="+t+", p="+p); xb = fx(r,t,p); yb = fy(r,t,p); zb = fz(r,t,p); err = Math.abs(x-xb)+Math.abs(y-yb)+Math.abs(z-zb); System.out.println("xb="+xb+", yb="+yb+", zb="+zb+", err="+err); rb = fr(xb,yb,zb); tb = ft(xb,yb,zb); pb = fp(xb,yb,zb); err = Math.abs(r-rb)+Math.abs(t-tb)+Math.abs(p-pb); System.out.println("rb="+rb+", tb="+tb+", pb="+pb+", err="+err); System.out.println(" "); System.out.println(" "); x = 1.5; y = -2.2; z = 3.1; System.out.println("x="+x+", y="+y+", z="+z); r = fr(x,y,z); t = ft(x,y,z); p = fp(x,y,z); System.out.println("r="+r+", t="+t+", p="+p); xb = fx(r,t,p); yb = fy(r,t,p); zb = fz(r,t,p); err = Math.abs(x-xb)+Math.abs(y-yb)+Math.abs(z-zb); System.out.println("xb="+xb+", yb="+yb+", zb="+zb+", err="+err); rb = fr(xb,yb,zb); tb = ft(xb,yb,zb); pb = fp(xb,yb,zb); err = Math.abs(r-rb)+Math.abs(t-tb)+Math.abs(p-pb); System.out.println("rb="+rb+", tb="+tb+", pb="+pb+", err="+err); System.out.println(" "); System.out.println(" "); x = 1.5; y = 2.2; z = -3.1; System.out.println("x="+x+", y="+y+", z="+z); r = fr(x,y,z); t = ft(x,y,z); p = fp(x,y,z); System.out.println("r="+r+", t="+t+", p="+p); xb = fx(r,t,p); yb = fy(r,t,p); zb = fz(r,t,p); err = Math.abs(x-xb)+Math.abs(y-yb)+Math.abs(z-zb); System.out.println("xb="+xb+", yb="+yb+", zb="+zb+", err="+err); rb = fr(xb,yb,zb); tb = ft(xb,yb,zb); pb = fp(xb,yb,zb); err = Math.abs(r-rb)+Math.abs(t-tb)+Math.abs(p-pb); System.out.println("rb="+rb+", tb="+tb+", pb="+pb+", err="+err); System.out.println(" "); System.out.println(" "); err = 0.0; r = r0; for(int i=0; i0.1) { System.out.println("error ="+err1); System.out.println("error at r="+r+", t="+t+", p="+p); System.out.println("error at xb="+xb+", yb="+yb+", zb="+zb); System.out.println("error at rb="+rb+", tb="+tb+", pb="+pb); System.out.println(" "); } else { err = err + err1; } p = p + dp; } tv[j] = t; // ? tt t = t + dt; } rv[i] = r; r = r + dr; } System.out.println("U[2][2][2]="+U[2][2][2]); System.out.println("r="+rv[2]+", t="+tv[2]+", p="+pv[2]); xb = fx(rv[2],tv[2],pv[2]); yb = fy(rv[2],tv[2],pv[2]); zb = fz(rv[2],tv[2],pv[2]); System.out.println("x="+xb+", y="+yb+", z="+zb); System.out.println(" "); System.out.println("U[3][3][3]="+U[3][3][3]); System.out.println("r="+rv[3]+", t="+tv[3]+", p="+pv[3]); xb = fx(rv[3],tv[3],pv[3]); yb = fy(rv[3],tv[3],pv[3]); zb = fz(rv[3],tv[3],pv[3]); System.out.println("x="+xb+", y="+yb+", z="+zb); System.out.println(" "); System.out.println("U[4][4][4]="+U[4][4][4]); System.out.println("r="+rv[4]+", t="+tv[4]+", p="+pv[4]); xb = fx(rv[4],tv[4],pv[4]); yb = fy(rv[4],tv[4],pv[4]); zb = fz(rv[4],tv[4],pv[4]); System.out.println("x="+xb+", y="+yb+", z="+zb); System.out.println(" "); System.out.println("sum of errors on "+(nr*nt*np)+" tests="+err); System.out.println(" "); System.out.println("spherical_deriv.java ends."); } double Uf(double r, double t, double p) { return r + r*t + r*Math.sin(t)*p; } double fx(double r, double t, double p) // x = r sin p cos t { return r * Math.sin(p) * Math.cos(t); } double fy(double r, double t, double p) // y = r sin p sin t { return r * Math.sin(p) * Math.sin(t); } double fz(double r, double t, double p) // z = r cos p { return r * Math.cos(p); } double fr(double x, double y, double z) // r = sqrt(x^2 + y^2 + z^2) { return Math.sqrt(x*x+y*y+z*z); } double ft(double x, double y, double z) // t = atan(y/x) { if(Math.abs(x)+Math.abs(y)<1.0e-10) return 0.0; // z axis return Math.atan2(y,x); } double fp(double x, double y, double z) // p = acos(z/r) { return Math.acos(z/Math.sqrt(x*x+y*y+z*z)); } public static void main (String[] args) { new spherical_deriv(); } }