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