Converting from coordinates in longitude latitude to heliocentric coordinate system with coordinates in arcsec

Script

double lon = Double.parseDouble(longitude);
double lat = Double.parseDouble(latitude);
double b0 = Double.parseDouble(b0String);
double sunr = Double.parseDouble(sunrString);

double arcsec2rad = 4.84813681e-6;

lon = lon * Math.PI/180.0;  //convert to radians
lat = lat * Math.PI/180.0;  //convert to radians

double cb0 = Math.cos(b0);
double clat = Math.cos(lat);
double clon = Math.cos(lon);
double sb0 = Math.sin(b0);
double slat = Math.sin(lat);
double slon = Math.sin(lon);

double crho = clat * cb0 * clon + slat*sb0;
if(crho > 1.0) {
  crho  = 1.0;
}
if(crho < -1.0) {
  crho = -1.0;
}
double rho = Math.acos(crho);

//theta is the positive angle CCW from north
double theta = -1.0 * (Math.atan2(cb0*clat*slon, slat-sb0*crho));

double ssunr = Math.sin(sunr*arcsec2rad);
double rho1 = Math.asin(ssunr*Math.sin(rho));
for (int i=0; i<5;i++){
  rho1 = Math.asin(ssunr*Math.sin(rho+rho1));
}

y = rho1*Math.cos(theta)/arcsec2rad;
x = -1.0*(rho1*Math.sin(theta)/arcsec2rad);