| 678 | } |
| 679 | |
| 680 | GMT_LOCAL void gmtplot_linearx_oblgrid (struct GMT_CTRL *GMT, struct PSL_CTRL *PSL, double w, double e, double s, double n, double dval) { |
| 681 | /* x gridlines in oblique coordinates for all but the Oblique Mercator projection [which already is oblique] */ |
| 682 | double *x = NULL, *lon = NULL, *lat = NULL, *lat_obl = NULL, tval, p_cap, s_cap; |
| 683 | unsigned int idup = 0, i, j, k, nx, np, nc1 = 0, nc2, npc, np1; |
| 684 | bool cap = false; |
| 685 | gmt_M_unused(PSL); gmt_M_unused(w); gmt_M_unused(e); |
| 686 | |
| 687 | /* Ideally we should determine the w/e/s/n of the oblique coordinates but here we will simply |
| 688 | * create oblique coordinates for the full 0/360/-90/90, convert to regular coordinates and |
| 689 | * then truncate points outside the actual w/e/s/n */ |
| 690 | |
| 691 | /* Do we have duplicate e and w boundaries ? */ |
| 692 | p_cap = fabs (GMT->current.setting.map_polar_cap[0]); |
| 693 | s_cap = -p_cap; |
| 694 | idup = (gmt_M_is_azimuthal(GMT)) ? 1 : 0; |
| 695 | cap = !doubleAlmostEqual (p_cap, 90.0); /* true if we have a polar cap specified */ |
| 696 | tval = (n - s) * GMT->current.setting.map_line_step / GMT->current.map.height; |
| 697 | |
| 698 | nx = gmtlib_linear_array (GMT, 0.0, TWO_PI, D2R * dval, D2R * GMT->current.map.frame.axis[GMT_X].phase, &x); |
| 699 | np = gmtlib_linear_array (GMT, -90.0, 90.0, tval, 0.0, &lat_obl); |
| 700 | np1 = nc2 = np - 1; /* Nominal number of points in path */ |
| 701 | lon = gmt_M_memory (GMT, NULL, np+2, double); /* Allow 2 more slots for possibly inserted cap-latitudes */ |
| 702 | lat = gmt_M_memory (GMT, NULL, np+2, double); |
| 703 | for (i = 0; i < nx - idup; i++) { /* For each oblique meridian to draw */ |
| 704 | /* Create lon,lat arrays of oblique coordinates for this meridian */ |
| 705 | for (k = j = 0; k < np; k++, j++) { |
| 706 | gmtlib_iobl (GMT, &lon[j], &lat[j], x[i], D2R * lat_obl[k]); /* Get regular coordinates of this point */ |
| 707 | lon[j] *= R2D; lat[j] *= R2D; /* Convert back to degrees */ |
| 708 | if (lat_obl[k] < s_cap && k < np1 && lat_obl[k+1] > s_cap) { /* Must insert S pole cap latitude point */ |
| 709 | j++; gmtlib_iobl (GMT, &lon[j], &lat[j], x[i], D2R * s_cap); |
| 710 | lon[j] *= R2D; lat[j] *= R2D; /* Back to degrees */ |
| 711 | nc1 = j; |
| 712 | } |
| 713 | else if (lat_obl[k] < p_cap && k < np1 && lat_obl[k+1] > p_cap) { /* Must insert N pole cap latitude point */ |
| 714 | j++; gmtlib_iobl (GMT, &lon[j], &lat[j], x[i], D2R * p_cap); |
| 715 | lon[j] *= R2D; lat[j] *= R2D; /* Back to degrees */ |
| 716 | nc2 = j; |
| 717 | } |
| 718 | } |
| 719 | if (cap) { /* Only plot the line between the two polar caps */ |
| 720 | npc = nc2 - nc1 + 1; /* Number of points along meridian bounded by caps */ |
| 721 | if ((GMT->current.plot.n = gmt_geo_to_xy_line (GMT, &lon[nc1], &lat[nc1], npc)) == 0) continue; |
| 722 | } |
| 723 | else { /* No polar cap in effect, plot entire meridian */ |
| 724 | if ((GMT->current.plot.n = gmt_geo_to_xy_line (GMT, lon, lat, j)) == 0) continue; |
| 725 | } |
| 726 | gmt_plot_line (GMT, GMT->current.plot.x, GMT->current.plot.y, GMT->current.plot.pen, GMT->current.plot.n, PSL_LINEAR); |
| 727 | } |
| 728 | if (nx) gmt_M_free (GMT, x); |
| 729 | if (cap) { /* Draw the polar cap(s) meridians with a separate lon spacing */ |
| 730 | nx = gmtlib_linear_array (GMT, 0.0, TWO_PI, D2R * GMT->current.setting.map_polar_cap[1], D2R * GMT->current.map.frame.axis[GMT_X].phase, &x); |
| 731 | for (i = 0; i < nx - idup; i++) { |
| 732 | for (k = j = 0; k < np; k++, j++) { |
| 733 | gmtlib_iobl (GMT, &lon[j], &lat[j], x[i], D2R * lat_obl[k]); /* Get regular coordinates of this point */ |
| 734 | lon[j] *= R2D; lat[j] *= R2D; /* Back to degrees */ |
| 735 | if (lat_obl[k] < s_cap && k < np1 && lat_obl[k+1] > s_cap) { /* Must insert S pole cap latitude point */ |
| 736 | j++; gmtlib_iobl (GMT, &lon[j], &lat[j], x[i], D2R * s_cap); |
| 737 | lon[j] *= R2D; lat[j] *= R2D; /* Back to degrees */ |
no test coverage detected