Calculate Sun data from date

Script

import java.text.SimpleDateFormat;

  static int JGREG= 15 + 31*(10+12*1582);

  public static double toJulian(int[] ymd) {
   int year=ymd[0];
   int month=ymd[1]; // jan=1, feb=2,...
   int day=ymd[2];
   int julianYear = year;
   if (year < 0) julianYear++;
   int julianMonth = month;
   if (month > 2) {
     julianMonth++;
   }
   else {
     julianYear--;
     julianMonth += 13;
   }

   double julian = (java.lang.Math.floor(365.25 * julianYear)
        + java.lang.Math.floor(30.6001*julianMonth) + day + 1720995.0);
   if (day + 31 * (month + 12 * year) >= JGREG) {
     // change over to Gregorian calendar
     int ja = (int)(0.01 * julianYear);
     julian += 2 - ja + (0.25 * ja);
   }
  return java.lang.Math.floor(julian);
 }



double radeg = 180.0/Math.PI;

// Julian date:


int[] dateIn =new int[3];
SimpleDateFormat sdf = new SimpleDateFormat("yyyy-MM-dd'T'HH:mm:ss");
Calendar cal = Calendar.getInstance();
cal.setTime(sdf.parse(dateTime));
dateIn[0]=cal.get(Calendar.YEAR);
dateIn[1]=cal.get(Calendar.MONTH)+1;
dateIn[2]=cal.get(Calendar.DAY_OF_MONTH);

// Julian Day number
double jd = toJulian(dateIn);

// Julian Centuries from 1900.0:
double t = (jd - 2415020.0)/36525.0;


//Geometric Mean Longitude (deg):
double  mnl = 279.69668 + 36000.76892*t + 0.0003025*Math.pow(t,2);
mnl = mnl % 360;

// Mean anomaly (deg):
double  mna = 358.47583 + 35999.04975*t -
        (0.000150*Math.pow(t,2)) - (0.0000033*Math.pow(t,3));
mna = mna % 360;

//Eccentricity of orbit:
double  e = 0.01675104 - (0.0000418*t) - (0.000000126*Math.pow(t,2));



// Sun's equation of center (deg):
double c = (1.919460 - (0.004789*t) - (0.000014*Math.pow(t,2)))*Math.sin(mna/radeg) + 
      (0.020094 - 0.000100*t)*Math.sin(2*mna/radeg) + (0.000293*Math.sin(3*mna/radeg));
// Sun's true geometric longitude (deg)
double true_long = (mnl + c) % 360.0;
double omega =  259.18 - 1934.142*t;
double lamda = true_long - 0.00569;
double lamda2 = lamda -.00479*Math.sin(omega/radeg);
double ob1 = 23.452294 - 0.0130125*t - 0.00000164*Math.pow(t,2)+ 0.000000503*Math.pow(t,3);

double k = 74.3646 + 1.395833 * t;
double diff = (lamda -k)/radeg;
double i = 7.25;
double x = Math.atan(Math.cos(lamda2/radeg)*-1*Math.tan(ob1/radeg))*radeg;
double y = Math.atan(Math.cos(diff)*-1*Math.tan(i/radeg))*radeg;
double pa = x + y;

position_angle = pa;

b0 = Math.asin(Math.sin(diff)*Math.sin(i/radeg))*radeg;

double ta = (mna +c) % 360;
double dist = 1.0000002 *(1.0 - Math.pow(e,2))/(1.0 +e*Math.cos(ta/radeg));
sun_radius = 959.63/dist;