MCPcopy Create free account
hub / github.com/GenericMappingTools/gmt / gmtlib_lonpath

Function gmtlib_lonpath

src/gmt_map.c:7860–7943  ·  view source on GitHub ↗

! . */

Source from the content-addressed store, hash-verified

7858
7859/*! . */
7860uint64_t gmtlib_lonpath (struct GMT_CTRL *GMT, double lon, double lat1, double lat2, double **x, double **y) {
7861 size_t n_alloc = 0;
7862 uint64_t n, k;
7863 int n_try, pos;
7864 bool keep_trying;
7865 double dlat, dlat0, *tlon = NULL, *tlat = NULL, x0, x1, y0, y1, d, min_gap;
7866
7867 if (GMT->current.map.meridian_straight == 2) { /* Special non-sampling for gmtselect/grdlandmask */
7868 gmt_M_malloc2 (GMT, tlon, tlat, 2U, NULL, double);
7869 tlon[0] = tlon[1] = lon;
7870 tlat[0] = lat1; tlat[1] = lat2;
7871 *x = tlon;
7872 *y = tlat;
7873 return (2ULL);
7874 }
7875
7876 if (GMT->current.map.meridian_straight) { /* Easy, just a straight line connect via quarter-points */
7877 gmt_M_malloc2 (GMT, tlon, tlat, 5, &n_alloc, double);
7878 tlon[0] = tlon[1] = tlon[2] = tlon[3] = tlon[4] = lon;
7879 dlat = lat2 - lat1;
7880 tlat[0] = lat1; tlat[1] = lat1 + 0.25 * dlat; tlat[2] = lat1 + 0.5 * dlat;
7881 tlat[3] = lat1 + 0.75 * dlat; tlat[4] = lat2;
7882 *x = tlon;
7883 *y = tlat;
7884 return (n = n_alloc);
7885 }
7886
7887 /* Must do general case */
7888 n = 0;
7889 min_gap = 0.1 * GMT->current.setting.map_line_step;
7890 if ((n_alloc = lrint (ceil (fabs (lat2 - lat1) / GMT->current.map.dlat))) == 0) return (0);
7891
7892 n_alloc++; /* So n_alloc is at least 2 */
7893 dlat0 = (lat2 - lat1) / n_alloc;
7894 pos = (dlat0 > 0.0);
7895
7896 k = n_alloc; n_alloc = 0;
7897 gmt_M_malloc2 (GMT, tlon, tlat, k, &n_alloc, double);
7898
7899 tlon[0] = lon;
7900 tlat[0] = lat1;
7901 gmt_geo_to_xy (GMT, tlon[0], tlat[0], &x0, &y0);
7902 while ((pos && (tlat[n] < lat2)) || (!pos && (tlat[n] > lat2))) {
7903 n++;
7904 if (n == n_alloc-1) {
7905 n_alloc += GMT_SMALL_CHUNK;
7906 tlon = gmt_M_memory (GMT, tlon, n_alloc, double);
7907 tlat = gmt_M_memory (GMT, tlat, n_alloc, double);
7908 }
7909 n_try = 0;
7910 keep_trying = true;
7911 dlat = dlat0;
7912 tlon[n] = lon;
7913 do {
7914 n_try++;
7915 tlat[n] = tlat[n-1] + dlat;
7916 if (gmt_M_y_is_lat (GMT, GMT_IN) && fabs (tlat[n]) > 90.0) tlat[n] = copysign (90.0, tlat[n]);
7917 gmt_geo_to_xy (GMT, tlon[n], tlat[n], &x1, &y1);

Callers 5

gmtplot_map_lonlineFunction · 0.85
gmt_geo_polarcap_segmentFunction · 0.85
gmt_graticule_pathFunction · 0.85
gmtlib_map_pathFunction · 0.85

Calls 2

gmt_geo_to_xyFunction · 0.85
hypotFunction · 0.85

Tested by

no test coverage detected