sunangle.c – this external uses the date and time from the seconds external to show the height and azimuth position of the sun. also shows how to handle an input list (same as A_GIMME)
/* sunangle.c - find the angle and azimuth of the sun at any latitude and longitude */ #include "m_pd.h" #include#include #define lotsOplaces 100 // generic, big-enough string length #define PI 3.141592653589793 typedef struct _sunangle // defines our object's internal variables for each instance in a patch { t_object x_ob; t_outlet *p_outletH; t_outlet *p_outletA; float year, month, date, hour, min, sec; float latitude; // first input float longitude; // second input } t_sunangle; t_class *sunangle_class; // global pointer to the object class - so max can reference the object // these are prototypes for the methods that are defined below void sunangle_bang(t_sunangle *x); void sunangle_seconds(t_sunangle *x, t_symbol *selector, int argcount, t_atom *argvec); void sunangle_latitude(t_sunangle *x, float f); void sunangle_longitude(t_sunangle *x, float f); void *sunangle_new(float lati, float longi); // prototype functions double calcJD(float, float, double); double calcTimeJulianCent(double); double calcEquationOfTime(double); double radToDeg(double); double degToRad(double); double calcGeomMeanLongSun(double); double calcGeomMeanAnomalySun(double); double calcSunDeclination(double); double calcObliquityCorrection(double); double calcMeanObliquityOfEcliptic(double); double calcSunApparentLong(double); double calcSunTrueLong(double); double calcGeomMeanLongSun(double); double calcSunEqOfCenter(double); double calcEccentricityEarthOrbit(double); void calcHA(float,float,float, float, float, float, float, float, float *, float *); //-------------------------------------------------------------------------- void sunangle_bang(t_sunangle *x) // x = reference to this instance of the object { float H, A; calcHA(x->latitude, x->longitude, x->year, x->month, x->date, x->hour, x->min, x->sec, &H, &A); outlet_float(x->p_outletH, H); outlet_float(x->p_outletA, A); } void sunangle_seconds(t_sunangle *x, t_symbol *selector, int argcount, t_atom *argvec) { if(argcount != 6) post("input to sunangle should be \"year month day hour minute second\""); else { if(argvec[0].a_type == A_FLOAT) x->year = argvec[0].a_w.w_float; else { post("the first item in list should be a float"); return; } if(argvec[1].a_type == A_FLOAT) x->month = argvec[1].a_w.w_float; else { post("the second item in list should be a float"); return; } if(argvec[2].a_type == A_FLOAT) x->date = argvec[2].a_w.w_float; else { post("the third item in list should be a float"); return; } if(argvec[3].a_type == A_FLOAT) x->hour = argvec[3].a_w.w_float; else { post("the fourth item in list should be a float"); return; } if(argvec[4].a_type == A_FLOAT) x->min = argvec[4].a_w.w_float; else { post("the fifth item in list should be a float"); return; } if(argvec[5].a_type == A_FLOAT) x->sec = argvec[5].a_w.w_float; else { post("the sixth item in list should be a float"); return; } } sunangle_bang(x); } void sunangle_latitude(t_sunangle *x, float f) { x->latitude = f; } // this gets called when something goes into inlet 2 void sunangle_longitude(t_sunangle *x, float f) { x->longitude = f; } //-------------------------------------------------------------------------- void *sunangle_new(float lati, float longi) // n = int argument typed into object box (A_DEFLONG) -- defaults to 0 if no args are typed { t_sunangle *x; // local variable (pointer to a t_sunangle data structure) x = (t_sunangle *)pd_new(sunangle_class); // create a new instance of this object // add inputs and outputs x->p_outletH = outlet_new(&x->x_ob, gensym("float")); x->p_outletA = outlet_new(&x->x_ob, gensym("float")); x->latitude = lati; // set initial (default) left operand value in the instance's data structure x->longitude = longi; // set initial (default) right operand value (n = variable passed to sunangle_new) return(x); // return a reference to the object instance } //-------------------------------------------------------------------------- void sunangle_setup(void) { sunangle_class = class_new(gensym("sunangle"), (t_newmethod)sunangle_new, 0, sizeof(t_sunangle), 0, A_DEFFLOAT, A_DEFFLOAT, 0); // class_new() loads our external into pd's memory so it can be used in a patch // sunangle_new = object creation method defined above class_addbang(sunangle_class, (t_method)sunangle_bang); class_addlist(sunangle_class, (t_method)sunangle_seconds); class_addmethod(sunangle_class, (t_method)sunangle_latitude, gensym("latitude"), A_FLOAT, 0); class_addmethod(sunangle_class, (t_method)sunangle_longitude, gensym("longitude"), A_FLOAT, 0); } // below this line is the actual sun position code //-------------------------------------------------------------------------- void calcHA(float lat, float lon, float year, float month, float date, float hour, float min, float sec, float *H, float *A){// sec is type double for calculations // declare calculated things double day, JD, T, minTime, time_offset, tst, sha, theta, cosPhi, exoatmElevation, refractionCorrection, te, azimuth; double zenith; double azDenom; int TZ = 0, DS = 0; // calculate astronomical times day = date+(hour+(min+sec/60.)/60.)/24.; JD = calcJD(year,month,day); T = calcTimeJulianCent(JD); minTime = calcEquationOfTime(T); time_offset = minTime-4*(-lon)+60*(TZ-DS); tst = hour*60+min+sec/60+time_offset; sha = tst/4-180; if (sha < -180) sha += 360.0; theta = calcSunDeclination(T); cosPhi = sin(degToRad(lat))*sin(degToRad(theta))+cos(degToRad(lat))* cos(degToRad(theta))*cos(degToRad(sha)); // exoatmospheric elevation angle exoatmElevation = 90 - radToDeg(acos(cosPhi)); // rudimentary refraction correction if (exoatmElevation > 85.0) { refractionCorrection = 0.0; } else { te = tan(degToRad(exoatmElevation)); if (exoatmElevation > 5.0) { refractionCorrection = 58.1 / te - 0.07 / (te*te*te) + 0.000086 / (te*te*te*te*te); } else if (exoatmElevation > -0.575) { refractionCorrection = 1735.0 + exoatmElevation * (-518.2 + exoatmElevation * (103.4 + exoatmElevation * (-12.79 + exoatmElevation * 0.711) ) ); } else { refractionCorrection = -20.774 / te; } refractionCorrection = refractionCorrection / 3600.0; } // refraction-corrected elevation angle *H = exoatmElevation + refractionCorrection; // work on azimuthal concerns if (cosPhi > 1.0) { cosPhi = 1.0; } else if (cosPhi < -1.0) { cosPhi = -1.0; } zenith = radToDeg(acos(cosPhi)); azDenom = cos(degToRad(lat))*sin(degToRad(zenith)); if (fabs(azDenom) > 0.001) { double azRad = (( sin(degToRad(lat)) * cos(degToRad(zenith)) ) - sin(degToRad(theta))) / azDenom; if (fabs(azRad) > 1.0) { if (azRad < 0) { azRad = -1.0; } else { azRad = 1.0; } } azimuth = 180.0 - radToDeg(acos(azRad)); if (sha > 0.0) { azimuth = -azimuth; } } else { if (lat > 0.0) { azimuth = 180.0; } else { azimuth = 0.0; } } if (azimuth < 0.0) { azimuth += 360.0; } *A = azimuth; } double calcJD(float year, float month, double day){ double JD, A, B; if(month<=2) { year = year-1; month = month+12; } A = floor((double)(year)/100.0); // note cast of year B = 2 - A + floor(A/4); JD = floor(365.25*(year + 4716.0)) + floor(30.6001*(month+1)) + day + B - 1524.5; return JD; } double calcTimeJulianCent(double jd){ double T; T = (jd - 2451545.0)/36525.0; return T; } double calcEquationOfTime(double t){ double minTime, epsilon, l0, e, m, y, sin2l0, sinm, cos2l0, sin4l0; double sin2m, Etime; epsilon = calcObliquityCorrection(t); l0 = calcGeomMeanLongSun(t); e = calcEccentricityEarthOrbit(t); m = calcGeomMeanAnomalySun(t); y = tan(degToRad(epsilon)/2.0); y = y*y; sin2l0 = sin(2.0 * degToRad(l0)); sinm = sin(degToRad(m)); cos2l0 = cos(2.0 * degToRad(l0)); sin4l0 = sin(4.0 * degToRad(l0)); sin2m = sin(2.0 * degToRad(m)); Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m; minTime = radToDeg(Etime)*4.0; return minTime; } double radToDeg(double angleRad){ return (180.0 * angleRad / PI); } double degToRad(double angleDeg){ return (PI * angleDeg / 180.0); } double calcGeomMeanLongSun(double t){ double L0 = 280.46646 + t * (36000.76983 + 0.0003032 * t); while(L0 > 360.0) { L0 -= 360.0; } while(L0 < 0.0) { L0 += 360.0; } return L0; // in degrees } double calcEccentricityEarthOrbit(double t){ double e = 0.016708634 - t * (0.000042037 + 0.0000001267 * t); return e; // unitless } double calcGeomMeanAnomalySun(double t){ double M = 357.52911 + t * (35999.05029 - 0.0001537 * t); return M; // in degrees } double calcObliquityCorrection(double t){ double e0 = calcMeanObliquityOfEcliptic(t); double omega = 125.04 - 1934.136 * t; double e = e0 + 0.00256 * cos(degToRad(omega)); return e; // in degrees } double calcSunDeclination(double t){ double e = calcObliquityCorrection(t); double lambda = calcSunApparentLong(t); double sint = sin(degToRad(e)) * sin(degToRad(lambda)); double theta = radToDeg(asin(sint)); return theta; // in degrees } double calcMeanObliquityOfEcliptic(double t){ double seconds = 21.448 - t*(46.8150 + t*(0.00059 - t*(0.001813))); double e0 = 23.0 + (26.0 + (seconds/60.0))/60.0; return e0; // in degrees } double calcSunApparentLong(double t){ double o = calcSunTrueLong(t); double omega = 125.04 - 1934.136 * t; double lambda = o - 0.00569 - 0.00478 * sin(degToRad(omega)); return lambda; // in degrees } double calcSunTrueLong(double t){ double l0 = calcGeomMeanLongSun(t); double c = calcSunEqOfCenter(t); double O = l0 + c; return O; // in degrees } double calcSunEqOfCenter(double t){ double m = calcGeomMeanAnomalySun(t); double mrad = degToRad(m); double sinm = sin(mrad); double sin2m = sin(mrad+mrad); double sin3m = sin(mrad+mrad+mrad); double C = sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) + sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289; return C; // in degrees }