// Copyright (c) CNES 2008 // // This software is part of CelestLab, a CNES toolbox for Scilab // // This software is governed by the CeCILL license under French law and // abiding by the rules of distribution of free software. You can use, // modify and/ or redistribute the software under the terms of the CeCILL // license as circulated by CEA, CNRS and INRIA at the following URL // 'http://www.cecill.info'. // ----------------------------------------------------------- //> Elevation from any location on Earth of a body //> (planets except Earth or Moon). //> - Elevation is the angle from the location between the body direction and its //> projection onto the horizontal plane (tangent to Earth surface). //> - The location is on the ground and is defined by its geodetic longitude //> and latitude. //> //> The available bodies are: Mercury, Venus, Mars, Jupiter, Saturn, //> Uranus, Neptune, Pluto, Moon. //> The colored areas correspond to the periods of time when the body //> is above the horizon. //> The contours levels give the elevation of the Sun //> (0 deg, -6 deg: civil twilight, -12 deg: nautical twilight, //> -18 deg: astronomical twilight) //> //> NB: //> - UT1 is here considered the same as UTC // // Author: A. Lamy // ----------------------------------------------------------- // Sun elevation // lon, lat: geodetic longitude and latitude (rad) (1x1 or 1xN) // t: UTC time (1x1 or 1xN) // body: "Sun" or "Moon" function [elev] = body_elev(body, lon, lat, t) N = max(size(lon,2), size(lat,2), size(t,2)); lon = lon .* ones(1:N); lat = lat .* ones(1:N); t = t .* ones(1:N); ra = CL_op_locTime(t, 'lon', lon, 'ra'); pos = CL_co_ell2car([ra; lat; %CL_eqRad * ones(lat)]); vert = CL_co_sph2car([ra; lat; ones(lat)]); if (body == "Sun") pos_body = CL_eph_sun(t); // Earth -> Sun elseif (body == "Moon") pos_body = CL_eph_moon(t); // Earth -> Moon else pos_body = CL_eph_planet(body, t); // Sun -> Body pos_body = CL_fr_convert("EOD", "ECI", t, pos_body); pos_body = pos_body + CL_eph_sun(t); // Earth -> body end elev = %pi/2 - CL_vectAngle(vert, pos_body-pos); endfunction // --------------------------------------------------------- // find (full) body name from beginning of name // => [] if does not exist // NB: Use global variables "Body_names" // --------------------------------------------------------- function [name] = getBodyName(str) name = []; str = stripblanks(str, %t); k = find(str == part(Body_names, 1:length(str))); if (length(k) == 1); name = Body_names(k); end endfunction // ----------------------------------------------------------- // Parameters input // ----------------------------------------------------------- Body_names = [ "Mercury", "Venus", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto", "Moon"]; lon = 0 * %CL_deg2rad; // geodetic longitude lat = 45 * %CL_deg2rad; // geodetic latitude name = "Venus"; // body name opt = 1; // night date_now = CL_dat_now("cal"); year = date_now(1); // year number (current year) desc_param = list(.. CL_defParam("Geodetic longitude of location", lon, units=['rad', 'deg'], valid='$x >= 0 & $x <= 360'),.. CL_defParam("Geodetic latitude of location", lat, units=['rad', 'deg'], valid='$x >= -90 & $x <= 90'),.. CL_defParam("Body name", name, typ='s', valid="getBodyName($x)<>[]"), .. CL_defParam("Time period (1=night, 2=any time)", opt, accv = [1,2]), .. CL_defParam("Year", year, valid="$x >= 1900 & $x <= 2100") .. ); [lon, lat, body, opt, year] = CL_inputParam(desc_param); body = getBodyName(body); // ----------------------------------------------------------- // computations // ----------------------------------------------------------- t0 = CL_dat_cal2cjd(year,1,1); tf = CL_dat_cal2cjd(year+1,1,1); yearlen = tf-t0; t = t0 : tf; s = (t - t0) / yearlen; utmin = 0; utmax = 24; ut = linspace(utmin, utmax, 120); [T, Ut] = ndgrid(t, ut); elev_sun = matrix(body_elev("Sun", lon, lat, T(:)'+Ut(:)'/24), size(T)); elev_body = matrix(body_elev(body, lon, lat, T(:)'+Ut(:)'/24), size(T)); // ----------------------------------------------------------- // plot // ----------------------------------------------------------- f = scf(); f.visible="on"; f.immediate_drawing="off"; nbcol = 50; f.color_map = 0.3 + 0.7*jetcolormap(nbcol); col_bk = addcolor([0,0,0]); col_g1 = addcolor([1,1,1] * 0.1); col_g2 = addcolor([1,1,1] * 0.5); a=gca(); CL_g_tag(a, 0); x = s; // abscissa for plot y = ut; // ordinate for plot z = elev_body * %CL_rad2deg; if (opt == 1) I = find(elev_sun * %CL_rad2deg >= 1); z(I) = -999; end zmin = 0; zmax = max(z); colorbar(zmin,zmax,colminmax=[1,nbcol],fmt="%.1f"); Sgrayplot(x, y, z, colminmax=[1,nbcol], ... zminmax=[zmin, zmax], colout=[0,0]); z = elev_sun * %CL_rad2deg; // Sun elevation: standard values for twilight levels = [0, -12]; sublevels = [-6, -18]; contour2d(x, y, z, sublevels, style=col_g2*ones(sublevels)); CL_g_tag(a,2); contour2d(x, y, z, levels, style=col_g1*ones(levels)); CL_g_tag(a,1); // plot adjustments h = CL_g_select(a, "Text",2); CL_g_delete(h); h = CL_g_select(a, "Text",1); if (h <> []) CL_g_set(h, "text", string(strtod(h.text))); h.font_foreground = col_g1; h.font_size = 3; h.font_style = 8; end h = CL_g_select(a, "Polyline",1); h.thickness = 2; h = CL_g_select(a, "Polyline",2); h.thickness = 1; // titles a.x_label.text = "Day / Month (" + string(year) + ")"; a.y_label.text = "Time (UT)"; str = body + " elevation (deg)"; a.title.text = sprintf("%s - lon = %.1f deg, lat = %.1f deg", str, lon*%CL_rad2deg, lat*%CL_rad2deg); // ticks m = [1:12, 1]; d = ones(m); y = [year * ones(1:12), year+1]; t_ticks = CL_dat_cal2cjd(y, m, d); labels = "1/" + string(m); x_ticks = (t_ticks - t0) / yearlen; a.x_ticks = tlist(['ticks', 'locations', 'labels'], x_ticks, labels); y_ticks = 0:2:24; a.y_ticks = tlist(['ticks', 'locations', 'labels'], y_ticks, string(y_ticks)); a.data_bounds = [0, 0; 1, 24]; a.tight_limits = "on"; CL_g_stdaxes(a, colg=col_g2); a.sub_ticks = [0,1]; // no substicks f.immediate_drawing="on"; f.visible="on";