cdm/src/main/java/ucar/unidata/geoloc/projection/sat/GEOSTransform.java
public double[] earthToSat(double geographic_lon, double geographic_lat) {

    geographic_lat = geographic_lat * DEG_TO_RAD;
    geographic_lon = geographic_lon * DEG_TO_RAD;

    double geocentric_lat = Math.atan(((r_pol * r_pol) / (r_eq * r_eq)) * Math.tan(geographic_lat));

    double r_earth = r_pol / Math.sqrt(1.0 - ((r_eq * r_eq - r_pol * r_pol) / (r_eq * r_eq)) * Math.cos(geocentric_lat) * Math.cos(geocentric_lat));

    double r_1 = h - r_earth * Math.cos(geocentric_lat) * Math.cos(geographic_lon - sub_lon);
    double r_2 = -r_earth * Math.cos(geocentric_lat) * Math.sin(geographic_lon - sub_lon);
    double r_3 = r_earth * Math.sin(geocentric_lat);

    if (r_1 > h) { // often two geoid intersect points, use the closer one.
      return new double[]{Double.NaN, Double.NaN};
    }

    double lamda_sat = Double.NaN;
    double theta_sat = Double.NaN;

    if (scan_geom.equals(GEOS)) { // GEOS (eg. SEVIRI, MSG)  CGMS 03, 4.4.3.2, Normalized Geostationary Projection
      lamda_sat = Math.atan(-r_2 / r_1);
      theta_sat = Math.asin(r_3 / Math.sqrt(r_1 * r_1 + r_2 * r_2 + r_3 * r_3));
    } else if (scan_geom.equals(GOES)) { // GOES (eg. GOES-R ABI)
      lamda_sat = Math.asin(-r_2 / Math.sqrt(r_1 * r_1 + r_2 * r_2 + r_3 * r_3));
      theta_sat = Math.atan(r_3 / r_1);
    }

    return new double[]{lamda_sat, theta_sat};
  }
public double[] earthToSat(double geographic_lon, double geographic_lat) {

    geographic_lat = geographic_lat * DEG_TO_RAD;
    geographic_lon = geographic_lon * DEG_TO_RAD;

    double geocentric_lat = Math.atan(((r_pol * r_pol) / (r_eq * r_eq)) * Math.tan(geographic_lat));

    double r_earth = r_pol / Math.sqrt(1.0 - ((r_eq * r_eq - r_pol * r_pol) / (r_eq * r_eq)) * Math.cos(geocentric_lat) * Math.cos(geocentric_lat));

    double r_1 = h - r_earth * Math.cos(geocentric_lat) * Math.cos(geographic_lon - sub_lon);
    double r_2 = -r_earth * Math.cos(geocentric_lat) * Math.sin(geographic_lon - sub_lon);
    double r_3 = r_earth * Math.sin(geocentric_lat);

    if (r_1 > h) { // often two geoid intersect points, use the closer one.
      return new double[]{Double.NaN, Double.NaN};
    }

    double lamda_sat = Double.NaN;
    double theta_sat = Double.NaN;

    if (scan_geom.equals(GEOS)) { // GEOS (eg. SEVIRI, MSG)  CGMS 03, 4.4.3.2, Normalized Geostationary Projection
      lamda_sat = Math.atan(-r_2 / r_1);
      theta_sat = Math.asin(r_3 / Math.sqrt(r_1 * r_1 + r_2 * r_2 + r_3 * r_3));
    } else if (scan_geom.equals(GOES)) { // GOES (eg. GOES-R ABI)
      lamda_sat = Math.asin(-r_2 / Math.sqrt(r_1 * r_1 + r_2 * r_2 + r_3 * r_3));
      theta_sat = Math.atan(r_3 / r_1);
    }

    return new double[]{lamda_sat, theta_sat};
  }
Transform geographic Earth coordinates to satellite view angle coordinate system
also known as the "intermediate" coordinate system in CGMS Normalized Geostationary Projection.

@param geographic_lon longitude, units: degrees
@param geographic_lat latitude, units: degrees
@return (lamda, theta) units: radian. This is the (x,y) or (East-West, North_South) view angle.