home *** CD-ROM | disk | FTP | other *** search
- /*--------------------------------------------------------------------
- * The GMT-system: @(#)gmt_map.c 2.56 09 Aug 1995
- *
- * Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
- * See README file for copying and redistribution conditions.
- *--------------------------------------------------------------------*/
-
- /*
- *
- * G M T _ M A P . C
- *
- *- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- * gmt_map.c contains code related to coordinate transformation
- *
- * Map Transformation Setup Routines
- * These routines initializes the selected map transformation
- * The names and main function are listed below
- * NB! Note that the transformation function does not check that they are
- * passed valid lon,lat numbers. I.e asking for log10 scaling using values
- * <= 0 results in problems.
- *
- * Map_projections include functions that will set up the transformation
- * between xy and latlon for several map projections.
- *
- * A few of the core coordinate transformation functions are based on similar
- * FORTRAN routines written by Pat Manley, Doug Shearer, and Bill Haxby, and
- * have been rewritten in C and subsequently streamlined. The rest is based
- * on P. Snyder, "Map Projections - a working manual", USGS Prof paper 1395.
- *
- * Transformations supported (both forward and inverse):
- * Linear x/y[/z] scaling
- * Polar (theta, r) scaling
- * Mercator
- * Stereographic
- * Albers Equal-Area Conic
- * Lambert Conformal Conic
- * Cassini Cylindrical
- * Oblique Mercator
- * TM Transverse Mercator
- * UTM Universal Transverse Mercator
- * Lambert Azimuthal Equal-Area
- * Mollweide Equal-Area
- * Hammer-Aitoff Equal-Area
- * Sinusoidal Equal-Area
- * Winkel Tripel
- * Orthographic
- * Azimuthal Equidistant
- * Robinson
- * Eckert IV
- * Cylindrical Equal-area (e.g., Peters, Gall, Behrmann)
- * Cylindrical Equidistant (Plate Carree)
- *
- * The ellipsoid used is selectable by editing the .gmtdefaults in your
- * home directory. If no such file, create one by running gmtdefaults.
- *
- * Usage: Initialize system by calling map_setup (separate module), and
- * then just use geo_to_xy() and xy_to_geo() functions.
- *
- * Author: Paul Wessel
- * Date: 27-JUL-1992
- * Version: v2.1
- *
- *
- * Functions include:
- *
- * azim_2_angle : Converts azimuth to angle on the map
- * break_through : Checks if we cross a map boundary
- * clip_to_map : Force polygon points to be inside map
- * compact_line : Remove redundant pen movements
- * geo_to_xy : Generic lon/lat to x/y
- * geo_to_xy_line : Same for polygons
- * geoz_to_xy : Generic 3-D lon/lat/z to x/y
- * get_origin : Find origin of projection based on pole and 2nd point
- * get_rotate_pole : Find rotation pole based on two points on great circle
- * grd_forward : Forward map-transform grid matrix from lon/lat to x/y
- * grd_inverse : Inversely transform grid matrix from x/y to lon/lat
- * grdproject_init : Initialize parameters for grid transformations
- * great_circle_dist : Returns great circle distance in degrees
- * hammer : Hammer-Aitoff equal area projection
- * ihammer : Inverse Hammer-Aitoff equal area projection
- * ilamb : Inverse Lambert conformal conic projection
- * ilambeq : Inverse Lambert azimuthal equal area projection
- * ilinearxy : Inverse linear projection
- * imerc : Inverse Mercator projection
- * imollweide : Inverse Mollweide projection
- * init_three_D : Initializes parameters needed for 3-D plots
- * ioblmrc : Inverse oblique Mercator projection, spherical only
- * iortho : Inverse orthographic projection
- * iplrs : Inverse Polar stereographic projection
- * itranslin : Inverse linear x projection
- * itranslog10 : Inverse log10 x projection
- * itranspowx : Inverse pow x projection
- * itranspowy : Inverse pow y projection
- * itm : Inverse TM projection
- * lamb : Lambert conformal conic projection
- * lambeq : Lambert azimuthal equal area projection
- * radial_crossing : Determine map crossing in the Lambert azimuthal equal area projection
- * latpath : Return path between 2 points of equal latitide
- * left_boundary : Return left boundary in x-inches
- * linearxy : Linear xy projection
- * lon_inside : Accounts for wrap-around in longitudes and checks for inside
- * lonpath : Return path between 2 points of equal longitude
- * map_crossing : Generic function finds crossings between line and map boundary
- * map_outside : Generic function determines if we're outside map boundary
- * map_path : Return latpat or lonpath
- * map_setup : Initialize map projection
- * merc : Mercator projection
- * mollweide : Mollweide projection
- * ellipse_crossing : Find map crossings in the Mollweide projection
- * move_to_rect : Move an outside point straight in to nearest edge
- * oblmrc : Spherical oblique Mercator projection
- * ortho : Orthographic projection
- * pen_status : Determines if pen is up or down
- * plrs : Polar stereographic projection
- * polar_outside : Determines if a point is outside polar projection region
- * pole_rotate_forward : Compute positions from oblique coordinates
- * pole_rotate_inverse : Compute oblique coordinates
- * project3D : Convert lon/lat/z to xx/yy/zz
- * radial_clip : Clip path outside radial region
- * radial_outside : Determine if point is outside radial region
- * radial_overlap : Determine overlap, always TRUE for his projection
- * rect_clip : Clip to rectangular region
- * rect_crossing : Find crossling between line and rect region
- * rect_outside : Determine if point is outside rect region
- * rect_outside2 : Determine if point is outside rect region (azimuthal proj only)
- * rect_overlap : Determine overlap between rect regions
- * right_boundary : Return x value of right map boundary
- * translin : Linear x projection
- * translog10 : Log10 x projection
- * transpowx : Linear pow x projection
- * transpowy : Linear pow y projection
- * tm : TM projection
- * xy_search : Find xy map boundary
- * two_d_to_three_d : Convert xyz to xy for entire array
- * vhammer : Initialize Hammer-Aitoff equal area projection
- * vlamb : Initialize Lambert conformal conic projection
- * vlambeq : Initialize Lambert azimuthal equal area projection
- * vmerc : Initialize Mercator projection
- * vmollweide : Initialize Mollweide projection
- * vortho : Initialize Orthographic projection
- * vplrs : Initialize Polar stereographic projection
- * vtm : Initialize TM projection
- * wesn_clip: Clip polygon to wesn boundaries
- * wesn_crossing : Find crossing between line and lon/lat rectangle
- * wesn_outside : Determine if a point is outside a lon/lat rectangle
- * wesn_overlap : Determine overlap between lon/lat rectangles
- * wesn_search : Search for extreme coordinates
- * wrap_around_check : Check if line wraps around due to Greenwich
- * x_to_xx : Generic linear x projection
- * xx_to_x : Generic inverse linear x projection
- * xy_to_geo : Generic inverse x/y to lon/lat projection
- * xyz_to_xy : Generic xyz to xy projection
- * y_to_yy : Generic linear y projection
- * yy_to_y : Generic inverse linear y projection
- * z_to_zz : Generic linear z projection
- * zz_to_z : Generic inverse linear z projection
- */
-
- #include "gmt.h"
-
- #define PI_2 (2.0 * M_PI)
-
- /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- *
- * S E C T I O N 1 : M A P - T R A N S F O R M A T I O N S
- *
- * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
-
- /*
- * map_setup sets up the transformations for the chosen map projection.
- * The user must pass:
- * w,e,s,n,parameters[0] - parameters[np-1] (np = number of parameters), and project:
- * w,e,s,n defines the area in degrees.
- * project == LINEAR, POLAR, MERCATOR, STEREO, LAMBERT, OBLIQUE_MERC,
- * TM, ORTHO, AZ_EQDIST, LAMB_AZ_EQ, WINKEL, ROBINSON, CASSINI, ALBERS. ECKERT, CYL_EQ, CYL_EQDIST
- * For LINEAR, we may have LINEAR, LOG10, or POW
- *
- * parameters[0] through parameters[np-1] mean different things to the various
- * projections, as explained below. (np also varies, of course)
- *
- * LINEAR projection:
- * parameters[0] is inch (or cm)/x_user_unit.
- * parameters[1] is inch (or cm)/y_user_unit. If 0, then yscale = xscale.
- * parameters[2] is pow for x^pow (if POW is on).
- * parameters[3] is pow for y^pow (if POW is on).
- *
- * POLAR (r,theta) projection:
- * parameters[0] is inch (or cm)/x_user_unit (radius)
- *
- * MERCATOR projection:
- * parameters[0] is in inch (or cm)/degree_longitude @ equator OR 1:xxxxx OR map-width
- *
- * STEREOgraphic projection:
- * parameters[0] is longitude of pole
- * parameters[1] is latitude of pole
- * parameters[2] is radius in inches (or cm) from pole to the latitude specified
- * by parameters[3] OR 1:xxxxx OR map-width.
- *
- * LAMBERT projection (Conic):
- * parameters[0] is first standard parallell
- * parameters[1] is second standard parallell
- * parameters[2] is scale in inch (or cm)/degree along parallells OR 1:xxxxx OR map-width
- *
- * OBLIQUE MERCATOR projection:
- * parameters[0] is longitude of origin
- * parameters[1] is latitude of origin
- * Definition a:
- * parameters[2] is longitude of second point along oblique equator
- * parameters[3] is latitude of second point along oblique equator
- * parameters[4] is scale in inch (or cm)/degree along oblique equator OR 1:xxxxx OR map-width
- * Definition b:
- * parameters[2] is azimuth along oblique equator at origin
- * parameters[3] is scale in inch (or cm)/degree along oblique equator OR 1:xxxxx OR map-width
- * Definition c:
- * parameters[2] is longitude of pole of projection
- * parameters[3] is latitude of pole of projection
- * parameters[4] is scale in inch (cm)/degree along oblique equator OR 1:xxxxx OR map-width
- *
- * TRANSVERSE MERCATOR (TM) projection
- * parameters[0] is central meridian
- * parameters[1] is scale in inch (cm)/degree along this meridian OR 1:xxxxx OR map-width
- *
- * UNIVERSAL TRANSVERSE MERCATOR (UTM) projection
- * parameters[0] is UTM zone (0-60, use negative for S hemisphere)
- * parameters[1] is scale in inch (cm)/degree along this meridian OR 1:xxxxx OR map-width
- *
- * LAMBERT AZIMUTHAL EQUAL AREA projection:
- * parameters[0] is longitude of origin
- * parameters[1] is latitude of origin
- * parameters[2] is radius in inches (or cm) from pole to the latitude specified
- * by parameters[3] OR 1:xxxxx OR map-width.
- *
- * ORTHOGRAPHIC AZIMUTHAL projection:
- * parameters[0] is longitude of origin
- * parameters[1] is latitude of origin
- * parameters[2] is radius in inches (or cm) from pole to the latitude specified
- * by parameters[3] OR 1:xxxxx OR map-width.
- *
- * AZIMUTHAL EQUIDISTANCE projection:
- * parameters[0] is longitude of origin
- * parameters[1] is latitude of origin
- * parameters[2] is radius in inches (or cm) from pole to the latitude specified
- * by parameters[3] OR 1:xxxxx OR map-width.
- *
- * MOLLWEIDE EQUAL-AREA projection
- * parameters[0] is longitude of origin
- * parametres[1] is in inch (or cm)/degree_longitude @ equator OR 1:xxxxx OR map-width
- *
- * HAMMER-AITOFF EQUAL-AREA projection
- * parameters[0] is longitude of origin
- * parametres[1] is in inch (or cm)/degree_longitude @ equator OR 1:xxxxx OR map-width
- *
- * SINUSOIDAL EQUAL-AREA projection
- * parameters[0] is longitude of origin
- * parameters[1] is in inch (or cm)/degree_longitude @ equator OR 1:xxxxx OR map-width
- *
- * WINKEL TRIPEL MODIFIED AZIMUTHAL projection
- * parameters[0] is longitude of origin
- * parametres[1] is in inch (or cm)/degree_longitude @ equator OR 1:xxxxx OR map-width
- *
- * ROBINSON PSEUDOCYLINDRICAL projection
- * parameters[0] is longitude of origin
- * parametres[1] is in inch (or cm)/degree_longitude @ equator OR 1:xxxxx OR map-width
- *
- * CASSINI projection
- * parameters[0] is longitude of origin
- * parameters[1] is latitude of origin
- * parametres[2] is scale in inch (cm)/degree along this meridian OR 1:xxxxx OR map-width
- *
- * ALBERS projection (Conic):
- * parameters[0] is first standard parallell
- * parameters[1] is second standard parallell
- * parameters[2] is scale in inch (or cm)/degree along parallells OR 1:xxxxx OR map-width
- *
- * ECKERT IV projection:
- * parameters[0] is longitude of origin
- * parameters[1] is scale in inch (or cm)/degree along parallells OR 1:xxxxx OR map-width
- *
- * CYLINDRICAL EQUAL-AREA projections (Behrmann, Gall, Peters):
- * parameters[0] is longitude of origin
- * parameters[1] is the standard parallell
- * parameters[2] is scale in inch (or cm)/degree along parallells OR 1:xxxxx OR map-width
- *
- * Pointers to the correct map transformation fuctions will be set up so that
- * there are no if tests to determine which routine to call. These pointers
- * are forward and inverse, and are called from geo_to_xy and xy_to_geo.
- *
- */
-
- int map_setup (west, east, south, north)
- double west, east, south, north; {
- int i, search;
- double f;
-
- if (!project_info.region) d_swap (south, east); /* Got w/s/e/n, make into w/e/s/n */
-
- if (west == east && south == north) {
- fprintf (stderr, "%s: GMT Fatal Error: No region selected - Aborts!\n", gmt_program);
- exit (-1);
- }
-
- if (east < west) east += 360.0;
- project_info.w = west; project_info.e = east; project_info.s = south; project_info.n = north;
- if (project_info.gave_map_width) project_info.units_pr_degree = FALSE;
-
- /* Set up ellipse parameters for the selected ellipsoid */
-
- frame_info.check_side = frame_info.horizontal = FALSE;
- EQ_RAD = gmtdefs.ellipse[gmtdefs.ellipsoid].eq_radius;
- f = gmtdefs.ellipse[gmtdefs.ellipsoid].flattening;
- ECC2 = 2 * f - f * f;
- ECC4 = ECC2 * ECC2;
- ECC6 = ECC2 * ECC4;
- ECC = d_sqrt (ECC2);
- M_PR_DEG = PI_2 * EQ_RAD / 360.0;
- if (!z_forward) z_forward = (PFI) translin;
- if (!z_inverse) z_inverse = (PFI) itranslin;
- gmtdefs.n_lon_nodes = gmtdefs.n_lat_nodes = 0;
-
- switch (project_info.projection) {
-
- case LINEAR: /* Linear transformations */
-
- search = map_init_linear ();
- break;
-
- case POLAR: /* Both lon/lat are actually theta, radius */
-
- search = map_init_polar ();
- break;
-
- case MERCATOR: /* Standard Mercator projection */
-
- search = map_init_merc ();
- break;
-
- case STEREO: /* Stereographic projection */
-
- search = map_init_stereo ();
- break;
-
- case LAMBERT: /* Lambert Conformal Conic */
-
- search = map_init_lambert();
- break;
-
- case OBLIQUE_MERC: /* Oblique Mercator */
-
- search = map_init_oblique ();
- break;
-
- case TM: /* Transverse Mercator */
-
- search = map_init_tm ();
- break;
-
- case UTM: /* Universal Transverse Mercator */
-
- search = map_init_utm ();
- break;
-
- case LAMB_AZ_EQ: /* Lambert Azimuthal Equal-Area */
-
- search = map_init_lambeq ();
- break;
-
- case ORTHO: /* Orthographic Projection */
-
- search = map_init_ortho ();
- break;
-
- case AZ_EQDIST: /* Azimuthal Equal-Distance Projection */
-
- search = map_init_azeqdist ();
- break;
-
- case MOLLWEIDE: /* Mollweide Equal-Area */
-
- search = map_init_mollweide ();
- break;
-
- case HAMMER: /* Hammer-Aitoff Equal-Area */
-
- search = map_init_hammer ();
- break;
-
- case WINKEL: /* Winkel Tripel */
-
- search = map_init_winkel ();
- break;
-
- case ECKERT: /* Eckert VI */
-
- search = map_init_eckert ();
- break;
-
- case CYL_EQ: /* Cylindrical Equal-Area */
-
- search = map_init_cyleq ();
- break;
-
- case CYL_EQDIST: /* Cylindrical Equaidistant */
-
- search = map_init_cyleqdist ();
- break;
-
- case ROBINSON: /* Robinson */
-
- search = map_init_robinson ();
- break;
-
- case SINUSOIDAL: /* Sinusoidal Equal-Area */
-
- search = map_init_sinusoidal ();
- break;
-
- case CASSINI: /* Cassini cylindrical */
-
- search = map_init_cassini ();
- break;
-
- case ALBERS: /* Albers Equal-Area Conic */
-
- search = map_init_albers ();
- break;
-
- default: /* No projection selected, die */
- fprintf (stderr, "%s: GMT Fatal Error: No projection selected - Aborts!\n", gmt_program);
- exit (-1);
- break;
- }
-
- gmt_map_width = fabs (project_info.xmax - project_info.xmin);
- gmt_map_height = fabs (project_info.ymax - project_info.ymin);
- gmt_half_map_size = 0.5 * gmt_map_width;
-
- if (!gmtdefs.n_lon_nodes) gmtdefs.n_lon_nodes = rint (gmt_map_width / gmtdefs.line_step);
- if (!gmtdefs.n_lat_nodes) gmtdefs.n_lat_nodes = rint (gmt_map_height / gmtdefs.line_step);
-
- gmtdefs.dlon = (project_info.e - project_info.w) / gmtdefs.n_lon_nodes;
- gmtdefs.dlat = (project_info.n - project_info.s) / gmtdefs.n_lat_nodes;
-
- if (search) {
- wesn_search ();
- gmtdefs.dlon = (project_info.e - project_info.w) / gmtdefs.n_lon_nodes;
- gmtdefs.dlat = (project_info.n - project_info.s) / gmtdefs.n_lat_nodes;
- }
-
- for (i = 0; i < 3; i++) {
- if (!project_info.xyz_pos[i]) { /* Max intervals negative */
- frame_info.anot_int[i] = -frame_info.anot_int[i];
- frame_info.frame_int[i] = -frame_info.frame_int[i];
- frame_info.grid_int[i] = -frame_info.grid_int[i];
- }
- }
-
- init_three_D ();
-
- /* Default for overlay plots is no shifting */
-
- if (!project_info.x_off_supplied && gmtdefs.overlay) gmtdefs.x_origin = 0.0;
- if (!project_info.y_off_supplied && gmtdefs.overlay) gmtdefs.y_origin = 0.0;
-
- return (0);
- }
-
- int init_three_D () {
- int i, easy;
- double tilt_angle, x, y, x0, x1, x2, y0, y1, y2, zmin, zmax;
-
- project_info.three_D = (z_project.view_azimuth != 180.0 || z_project.view_elevation != 90.0);
- if (project_info.three_D) gmtdefs.basemap_type = 1; /* Only plain frame for 3-D */
-
- project_info.z_scale = project_info.z_pars[0];
- if (project_info.z_scale < 0.0) project_info.xyz_pos[2] = FALSE; /* User wants z to increase down */
- if (project_info.z_level == MAXDOUBLE) project_info.z_level = (project_info.xyz_pos[2]) ? project_info.z_bottom : project_info.z_top;
-
- switch (project_info.xyz_projection[2]) {
- case LINEAR: /* Regular scaling */
- zmin = (project_info.xyz_pos[2]) ? project_info.z_bottom : project_info.z_top;
- zmax = (project_info.xyz_pos[2]) ? project_info.z_top : project_info.z_bottom;
- z_forward = (PFI) translin;
- z_inverse = (PFI) itranslin;
- break;
- case LOG10: /* Log10 transformation */
- if (project_info.z_bottom <= 0.0 || project_info.z_top <= 0.0) {
- fprintf (stderr, "%s: GMT SYNTAX ERROR for -Jz -JZ option: limits must be positive for log10 projection\n", gmt_program);
- exit (-1);
- }
- zmin = (project_info.xyz_pos[2]) ? d_log10 (project_info.z_bottom) : d_log10 (project_info.z_top);
- zmax = (project_info.xyz_pos[2]) ? d_log10 (project_info.z_top) : d_log10 (project_info.z_bottom);
- z_forward = (PFI) translog10;
- z_inverse = (PFI) itranslog10;
- break;
- case POW: /* x^y transformation */
- project_info.xyz_pow[2] = project_info.z_pars[1];
- zmin = (project_info.xyz_pos[2]) ? pow (project_info.z_bottom, project_info.xyz_pow[2]) : pow (project_info.z_top, project_info.xyz_pow[2]);
- zmax = (project_info.xyz_pos[2]) ? pow (project_info.z_top, project_info.xyz_pow[2]) : pow (project_info.z_bottom, project_info.xyz_pow[2]);
- z_forward = (PFI) transpowz;
- z_inverse = (PFI) itranspowz;
- }
- if (project_info.compute_scale[2]) project_info.z_scale /= fabs (zmin - zmax);
- project_info.zmax = (zmax - zmin) * project_info.z_scale;
- project_info.z0 = -zmin * project_info.z_scale;
-
- if (z_project.view_azimuth >= 360.0) z_project.view_azimuth -= 360.0;
- z_project.quadrant = ceil (z_project.view_azimuth / 90.0);
- z_project.view_azimuth -= 180.0; /* Turn into direction instead */
- if (z_project.view_azimuth < 0.0) z_project.view_azimuth += 360.0;
- z_project.view_azimuth *= D2R;
- z_project.view_elevation *= D2R;
- z_project.cos_az = cos (z_project.view_azimuth);
- z_project.sin_az = sin (z_project.view_azimuth);
- z_project.cos_el = cos (z_project.view_elevation);
- z_project.sin_el = sin (z_project.view_elevation);
- geoz_to_xy (project_info.w, project_info.s, project_info.z_bottom, &x0, &y0);
- geoz_to_xy (project_info.e, project_info.s, project_info.z_bottom, &x1, &y1);
- geoz_to_xy (project_info.w, project_info.n, project_info.z_bottom, &x2, &y2);
- z_project.phi[0] = d_atan2 (y1 - y0, x1 - x0) * R2D;
- z_project.phi[1] = d_atan2 (y2 - y0, x2 - x0) * R2D;
- z_project.phi[2] = 90.0;
- tilt_angle = (z_project.phi[0] + 90.0 - z_project.phi[1]) * D2R;
- z_project.k = (fabs (z_project.cos_az) > fabs (z_project.sin_az)) ? 0 : 1;
- z_project.xshrink[0] = hypot (z_project.cos_az, z_project.sin_az * z_project.sin_el);
- z_project.yshrink[0] = hypot (z_project.sin_az, z_project.cos_az * z_project.sin_el);
- z_project.xshrink[1] = z_project.yshrink[0];
- z_project.yshrink[1] = z_project.xshrink[0];
- z_project.yshrink[0] *= fabs (cos (tilt_angle)); /* New */
- z_project.yshrink[1] *= fabs (cos (tilt_angle));
- z_project.xshrink[2] = fabs (z_project.cos_el);
- z_project.yshrink[2] = (fabs (z_project.cos_az) > fabs (z_project.sin_az)) ? fabs (z_project.cos_az) : fabs (z_project.sin_az);
- z_project.tilt[0] = tan (tilt_angle);
- z_project.tilt[1] = -z_project.tilt[0];
- z_project.tilt[2] = (fabs (z_project.cos_az) > fabs (z_project.sin_az)) ? tan (-z_project.phi[0] * D2R) : tan (-z_project.phi[1] * D2R);
-
- /* Determine min/max y of plot */
-
- /* easy = TRUE means we can use 4 corner points to find min x and y, FALSE
- measn we must search along wesn perimeter the hard way */
-
- switch (project_info.projection) {
- case LINEAR:
- case POLAR:
- case MERCATOR:
- case OBLIQUE_MERC:
- case CYL_EQ:
- case CYL_EQDIST:
- easy = TRUE;
- break;
- case LAMBERT:
- case TM:
- case UTM:
- case CASSINI:
- case STEREO:
- case ALBERS:
- case LAMB_AZ_EQ:
- case ORTHO:
- case AZ_EQDIST:
- case SINUSOIDAL:
- case MOLLWEIDE:
- case HAMMER:
- case WINKEL:
- case ECKERT:
- case ROBINSON:
- easy = (!project_info.region);
- break;
- default:
- easy = FALSE;
- break;
- }
-
- if (!project_info.three_D) easy = TRUE;
-
- z_project.corner_x[0] = z_project.corner_x[3] = (project_info.xyz_pos[0]) ? project_info.w : project_info.e;
- z_project.corner_x[1] = z_project.corner_x[2] = (project_info.xyz_pos[0]) ? project_info.e : project_info.w;
- z_project.corner_y[0] = z_project.corner_y[1] = (project_info.xyz_pos[1]) ? project_info.s : project_info.n;
- z_project.corner_y[2] = z_project.corner_y[3] = (project_info.xyz_pos[1]) ? project_info.n : project_info.s;
- z_project.xmin = z_project.ymin = MAXDOUBLE; z_project.xmax = z_project.ymax = -MAXDOUBLE;
-
- if (easy) {
- double xx[4], yy[4];
-
- xx[0] = xx[3] = project_info.xmin; xx[1] = xx[2] = project_info.xmax;
- yy[0] = yy[1] = project_info.ymin; yy[2] = yy[3] = project_info.ymax;
-
- for (i = 0; i < 4; i++) {
- xy_do_z_to_xy (xx[i], yy[i], project_info.z_bottom, &x, &y);
- z_project.ymin = MIN (z_project.ymin, y);
- z_project.ymax = MAX (z_project.ymax, y);
- z_project.xmin = MIN (z_project.xmin, x);
- z_project.xmax = MAX (z_project.xmax, x);
- xy_do_z_to_xy (xx[i], yy[i], project_info.z_top, &x, &y);
- z_project.ymin = MIN (z_project.ymin, y);
- z_project.ymax = MAX (z_project.ymax, y);
- z_project.xmin = MIN (z_project.xmin, x);
- z_project.xmax = MAX (z_project.xmax, x);
- }
- }
- else {
- for (i = 0; i < gmtdefs.n_lon_nodes; i++) { /* S and N */
- geoz_to_xy (project_info.w + i * gmtdefs.dlon, project_info.s, project_info.z_bottom, &x, &y);
- z_project.ymin = MIN (z_project.ymin, y);
- z_project.ymax = MAX (z_project.ymax, y);
- z_project.xmin = MIN (z_project.xmin, x);
- z_project.xmax = MAX (z_project.xmax, x);
- if (project_info.z_top != project_info.z_bottom) {
- geoz_to_xy (project_info.w + i * gmtdefs.dlon, project_info.s, project_info.z_top, &x, &y);
- z_project.ymin = MIN (z_project.ymin, y);
- z_project.ymax = MAX (z_project.ymax, y);
- z_project.xmin = MIN (z_project.xmin, x);
- z_project.xmax = MAX (z_project.xmax, x);
- }
- geoz_to_xy (project_info.w + i * gmtdefs.dlon, project_info.n, project_info.z_bottom, &x, &y);
- z_project.ymin = MIN (z_project.ymin, y);
- z_project.ymax = MAX (z_project.ymax, y);
- z_project.xmin = MIN (z_project.xmin, x);
- z_project.xmax = MAX (z_project.xmax, x);
- if (project_info.z_top != project_info.z_bottom) {
- geoz_to_xy (project_info.w + i * gmtdefs.dlon, project_info.n, project_info.z_top, &x, &y);
- z_project.ymin = MIN (z_project.ymin, y);
- z_project.ymax = MAX (z_project.ymax, y);
- z_project.xmin = MIN (z_project.xmin, x);
- z_project.xmax = MAX (z_project.xmax, x);
- }
- }
- for (i = 0; i < gmtdefs.n_lat_nodes; i++) { /* W and E */
- geoz_to_xy (project_info.w, project_info.s + i * gmtdefs.dlat, project_info.z_bottom, &x, &y);
- z_project.ymin = MIN (z_project.ymin, y);
- z_project.ymax = MAX (z_project.ymax, y);
- z_project.xmin = MIN (z_project.xmin, x);
- z_project.xmax = MAX (z_project.xmax, x);
- if (project_info.z_top != project_info.z_bottom) {
- geoz_to_xy (project_info.w, project_info.s + i * gmtdefs.dlat, project_info.z_top, &x, &y);
- z_project.ymin = MIN (z_project.ymin, y);
- z_project.ymax = MAX (z_project.ymax, y);
- z_project.xmin = MIN (z_project.xmin, x);
- z_project.xmax = MAX (z_project.xmax, x);
- }
- geoz_to_xy (project_info.e, project_info.s + i * gmtdefs.dlat, project_info.z_bottom, &x, &y);
- z_project.ymin = MIN (z_project.ymin, y);
- z_project.ymax = MAX (z_project.ymax, y);
- z_project.xmin = MIN (z_project.xmin, x);
- z_project.xmax = MAX (z_project.xmax, x);
- if (project_info.z_top != project_info.z_bottom) {
- geoz_to_xy (project_info.e, project_info.s + i * gmtdefs.dlat, project_info.z_top, &x, &y);
- z_project.ymin = MIN (z_project.ymin, y);
- z_project.ymax = MAX (z_project.ymax, y);
- z_project.xmin = MIN (z_project.xmin, x);
- z_project.xmax = MAX (z_project.xmax, x);
- }
- }
- }
-
- z_project.face[0] = (z_project.quadrant == 1 || z_project.quadrant == 2) ? 0 : 1;
- z_project.face[1] = (z_project.quadrant == 1 || z_project.quadrant == 4) ? 2 : 3;
- z_project.face[2] = (z_project.view_elevation >= 0.0) ? 4 : 5;
- z_project.draw[0] = (z_project.quadrant == 1 || z_project.quadrant == 4) ? TRUE : FALSE;
- z_project.draw[1] = (z_project.quadrant == 3 || z_project.quadrant == 4) ? TRUE : FALSE;
- z_project.draw[2] = (z_project.quadrant == 2 || z_project.quadrant == 3) ? TRUE : FALSE;
- z_project.draw[3] = (z_project.quadrant == 1 || z_project.quadrant == 2) ? TRUE : FALSE;
- z_project.sign[0] = z_project.sign[3] = -1.0;
- z_project.sign[1] = z_project.sign[2] = 1.0;
- z_project.z_axis = (z_project.quadrant%2) ? z_project.quadrant : z_project.quadrant - 2;
- }
-
- /*
- * GENERIC TRANSFORMATION ROUTINES FOR THE LINEAR PROJECTION
- */
-
- int x_to_xx (x, xx)
- double x, *xx; {
- /* Converts x to xx using the current linear projection */
-
- (*x_forward) (x, xx);
- (*xx) = (*xx) * project_info.x_scale + project_info.x0;
- }
-
- int y_to_yy (y, yy)
- double y, *yy; {
- /* Converts y to yy using the current linear projection */
-
- (*y_forward) (y, yy);
- (*yy) = (*yy) * project_info.y_scale + project_info.y0;
- }
-
- int z_to_zz (z, zz)
- double z, *zz; {
- /* Converts z to zz using the current linear projection */
-
- (*z_forward) (z, zz);
- (*zz) = (*zz) * project_info.z_scale + project_info.z0;
- }
-
- int xx_to_x (x, xx)
- double *x, xx; {
- /* Converts xx back to x using the current linear projection */
-
- xx = (xx - project_info.x0) / project_info.x_scale;
-
- (*x_inverse) (x, xx);
- }
-
- int yy_to_y (y, yy)
- double *y, yy; {
- /* Converts yy back to y using the current linear projection */
-
- yy = (yy - project_info.y0) / project_info.y_scale;
-
- (*y_inverse) (y, yy);
- }
-
- int zz_to_z (z, zz)
- double *z, zz; {
- /* Converts zz back to z using the current linear projection */
-
- zz = (zz - project_info.z0) / project_info.z_scale;
-
- (*z_inverse) (z, zz);
- }
-
- int geo_to_xy (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Converts lon/lat to x/y using the current projection */
-
- (*forward) (lon, lat, x, y);
- (*x) = (*x) * project_info.x_scale + project_info.x0;
- (*y) = (*y) * project_info.y_scale + project_info.y0;
- }
-
- int xy_to_geo (lon, lat, x, y)
- double *lon, *lat, x, y; {
- /* Converts x/y to lon/lat using the current projection */
-
- x = (x - project_info.x0) / project_info.x_scale;
- y = (y - project_info.y0) / project_info.y_scale;
-
- (*inverse) (lon, lat, x, y);
- }
-
- int geoz_to_xy (x, y, z, x_out, y_out)
- double x, y, z;
- double *x_out, *y_out; { /* Map-projects xy first, the projects xyz onto xy plane */
- double x0, y0, z0;
- geo_to_xy (x, y, &x0, &y0);
- z_to_zz (z, &z0);
- *x_out = x0 * z_project.cos_az - y0 * z_project.sin_az;
- *y_out = (x0 * z_project.sin_az + y0 * z_project.cos_az) * z_project.sin_el + z0 * z_project.cos_el;
- }
-
- int xyz_to_xy (x, y, z, x_out, y_out)
- double x, y, z;
- double *x_out, *y_out; { /* projects xyz onto xy plane */
- *x_out = x * z_project.cos_az - y * z_project.sin_az;
- *y_out = (x * z_project.sin_az + y * z_project.cos_az) * z_project.sin_el + z * z_project.cos_el;
- }
-
- int xy_do_z_to_xy (x, y, z, x_out, y_out)
- double x, y, z;
- double *x_out, *y_out; { /* projects xy (inches) z (z-value) onto xy plane */
- double z_out;
-
- z_to_zz (z, &z_out);
- *x_out = x * z_project.cos_az - y * z_project.sin_az;
- *y_out = (x * z_project.sin_az + y * z_project.cos_az) * z_project.sin_el + z_out * z_project.cos_el;
- }
-
- int project3D (x, y, z, x_out, y_out, z_out)
- double x, y, z, *x_out, *y_out, *z_out; {
- geo_to_xy (x, y, x_out, y_out);
- z_to_zz (z, z_out);
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE LINEAR PROJECTION
- */
-
- int map_init_linear () {
- BOOLEAN degree;
- double xmin, xmax, ymin, ymax;
-
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- forward = (PFI) linearxy;
- inverse = (PFI) ilinearxy;
- degree = (project_info.pars[4] == 1.0);
- if (degree) {
- project_info.central_meridian = 0.5 * (project_info.w + project_info.e);
- gmt_world_map = (fabs (project_info.w - project_info.e) == 360.0);
- }
- project_info.x_scale = project_info.pars[0];
- project_info.y_scale = project_info.pars[1];
- if (project_info.x_scale < 0.0) project_info.xyz_pos[0] = FALSE; /* User wants x to increase left */
- if (project_info.y_scale < 0.0) project_info.xyz_pos[1] = FALSE; /* User wants y to increase down */
-
- switch (project_info.xyz_projection[0]) {
- case LINEAR: /* Regular scaling */
- x_forward = (PFI) ((degree) ? translind : translin);
- x_inverse = (PFI) itranslin;
- if (project_info.xyz_pos[0]) {
- (*x_forward) (project_info.w, &xmin);
- (*x_forward) (project_info.e, &xmax);
- }
- else {
- (*x_forward) (project_info.e, &xmin);
- (*x_forward) (project_info.w, &xmax);
- }
- break;
- case LOG10: /* Log10 transformation */
- if (project_info.w <= 0.0 || project_info.e <= 0.0) {
- fprintf (stderr, "%s: GMT SYNTAX ERROR -Jx option: Limits must be positive for log10 option\n", gmt_program);
- exit (-1);
- }
- xmin = (project_info.xyz_pos[0]) ? d_log10 (project_info.w) : d_log10 (project_info.e);
- xmax = (project_info.xyz_pos[0]) ? d_log10 (project_info.e) : d_log10 (project_info.w);
- x_forward = (PFI) translog10;
- x_inverse = (PFI) itranslog10;
- break;
- case POW: /* x^y transformation */
- project_info.xyz_pow[0] = project_info.pars[2];
- project_info.xyz_ipow[0] = 1.0 / project_info.pars[2];
- xmin = (project_info.xyz_pos[0]) ? pow (project_info.w, project_info.xyz_pow[0]) : pow (project_info.e, project_info.xyz_pow[0]);
- xmax = (project_info.xyz_pos[0]) ? pow (project_info.e, project_info.xyz_pow[0]) : pow (project_info.w, project_info.xyz_pow[0]);
- x_forward = (PFI) transpowx;
- x_inverse = (PFI) itranspowx;
- break;
- }
- switch (project_info.xyz_projection[1]) {
- case LINEAR: /* Regular scaling */
- ymin = (project_info.xyz_pos[1]) ? project_info.s : project_info.n;
- ymax = (project_info.xyz_pos[1]) ? project_info.n : project_info.s;
- y_forward = (PFI) translin;
- y_inverse = (PFI) itranslin;
- break;
- case LOG10: /* Log10 transformation */
- if (project_info.s <= 0.0 || project_info.n <= 0.0) {
- fprintf (stderr, "%s: GMT SYNTAX ERROR -Jx option: Limits must be positive for log10 option\n", gmt_program);
- exit (-1);
- }
- ymin = (project_info.xyz_pos[1]) ? d_log10 (project_info.s) : d_log10 (project_info.n);
- ymax = (project_info.xyz_pos[1]) ? d_log10 (project_info.n) : d_log10 (project_info.s);
- y_forward = (PFI) translog10;
- y_inverse = (PFI) itranslog10;
- break;
- case POW: /* x^y transformation */
- project_info.xyz_pow[1] = project_info.pars[3];
- project_info.xyz_ipow[1] = 1.0 / project_info.pars[3];
- ymin = (project_info.xyz_pos[1]) ? pow (project_info.s, project_info.xyz_pow[1]) : pow (project_info.n, project_info.xyz_pow[1]);
- ymax = (project_info.xyz_pos[1]) ? pow (project_info.n, project_info.xyz_pow[1]) : pow (project_info.s, project_info.xyz_pow[1]);
- y_forward = (PFI) transpowy;
- y_inverse = (PFI) itranspowy;
- }
-
- /* Was given axes length instead of scale? */
-
- if (project_info.compute_scale[0]) project_info.x_scale /= fabs (xmin - xmax);
- if (project_info.compute_scale[1]) project_info.y_scale /= fabs (ymin - ymax);
-
- map_setxy (xmin, xmax, ymin, ymax);
- outside = (PFI) rect_outside;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- frame_info.horizontal = TRUE;
-
- return (FALSE);
- }
-
- int linearxy (x, y, x_i, y_i)
- double x, y, *x_i, *y_i; { /* Transform both x and y linearly */
- (*x_forward) (x, x_i);
- (*y_forward) (y, y_i);
- }
-
- int ilinearxy (x, y, x_i, y_i)
- double *x, *y, x_i, y_i; { /* Inversely transform both x and y linearly */
- (*x_inverse) (x, x_i);
- (*y_inverse) (y, y_i);
- }
-
- int translin (forw, inv) /* Linear forward */
- double forw, *inv; {
- *inv = forw;
- }
-
- int translind (forw, inv) /* Linear forward, but with degrees*/
- double forw, *inv; {
- while ((forw - project_info.central_meridian) < -180.0) forw += 360.0;
- while ((forw - project_info.central_meridian) > 180.0) forw -= 360.0;
- *inv = forw - project_info.central_meridian;
- }
-
- int itranslin (forw, inv) /* Linear inverse */
- double *forw, inv; {
- *forw = inv;
- }
-
- int translog10 (forw, inv) /* Log10 forward */
- double forw, *inv; {
- *inv = d_log10 (forw);
- }
-
- int itranslog10 (forw, inv) /* Log10 inverse */
- double *forw, inv; {
- *forw = pow (10.0, inv);
- }
-
- int transpowx (x, x_in) /* pow x forward */
- double x, *x_in; {
- *x_in = pow (x, project_info.xyz_pow[0]);
- }
-
- int itranspowx (x, x_in) /* pow x inverse */
- double *x, x_in; {
- *x = pow (x_in, project_info.xyz_ipow[0]);
- }
-
- int transpowy (y, y_in) /* pow y forward */
- double y, *y_in; {
- *y_in = pow (y, project_info.xyz_pow[1]);
- }
-
- int itranspowy (y, y_in) /* pow y inverse */
- double *y, y_in; {
- *y = pow (y_in, project_info.xyz_ipow[1]);
- }
-
- int transpowz (z, z_in) /* pow z forward */
- double z, *z_in; {
- *z_in = pow (z, project_info.xyz_pow[2]);
- }
-
- int itranspowz (z, z_in) /* pow z inverse */
- double *z, z_in; {
- *z = pow (z_in, project_info.xyz_ipow[2]);
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR POLAR (theta,r) PROJECTION
- */
-
- int map_init_polar () {
- double xmin, xmax, ymin, ymax;
-
- vpolar (0.5 * (project_info.w + project_info.e));
- /* project_info.w = 0.0; project_info.e = 360.0; */
- /* if (project_info.gave_map_width) project_info.pars[0] *= 0.5 / project_info.n; */ /* north is max radius */
- if (project_info.s == 0.0) project_info.edge[0] = FALSE;
- if (fabs (project_info.e - project_info.w) == 360.0) project_info.edge[1] = project_info.edge[3] = FALSE;
- left_edge = (PFD) left_circle;
- right_edge = (PFD) right_circle;
- forward = (PFI) polar;
- inverse = (PFI) ipolar;
- gmt_world_map = TRUE;
- xy_search (&xmin, &xmax, &ymin, &ymax);
- project_info.x_scale = project_info.y_scale = project_info.pars[0];
- map_setinfo (xmin, xmax, ymin, ymax, project_info.pars[0]);
- /* xmin = ymin = -project_info.n; xmax = ymax = project_info.n;
- project_info.x_scale = project_info.y_scale = project_info.pars[0]; */
-
- /* project_info.r = 0.5 * project_info.xmax; */
- project_info.r = project_info.y_scale * project_info.n;
- outside = (PFI) polar_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- frame_info.horizontal = TRUE;
- gmtdefs.degree_format = 6; /* Special labeling case */
-
- return (FALSE);
- }
-
- int vpolar (lon0)
- double lon0; {
- /* Set up a Polar (theta,r) transformation */
-
- project_info.central_meridian = lon0;
- project_info.north_pole = FALSE;
- return (0);
- }
-
- int polar (x, y, x_i, y_i)
- double x, y, *x_i, *y_i; { /* Transform both x and y linearly */
- *x_i = y * cosd (x);
- *y_i = y * sind (x);
- }
-
- int ipolar (x, y, x_i, y_i)
- double *x, *y, x_i, y_i; { /* Inversely transform both x and y linearly */
- *x = R2D * d_atan2 (y_i, x_i);
- *y = hypot (x_i, y_i);
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE MERCATOR PROJECTION
- */
-
- int map_init_merc () {
- double xmin, xmax, ymin, ymax;
-
- if (project_info.s == -90.0 || project_info.n == 90.0) {
- fprintf (stderr, "%s: GMT SYNTAX ERROR -R option: Cannot include south/north poles with Mercator projection!\n", gmt_program);
- exit (-1);
- }
-
- vmerc (0.5 * (project_info.w + project_info.e));
- merc (project_info.w, project_info.s, &xmin, &ymin);
- merc (project_info.e, project_info.n, &xmax, &ymax);
- if (project_info.gave_map_width) project_info.pars[0] /= (M_PR_DEG * fabs (project_info.w - project_info.e));
- if (project_info.units_pr_degree) project_info.pars[0] /= M_PR_DEG;
- project_info.x_scale = project_info.y_scale = project_info.pars[0];
- map_setxy (xmin, xmax, ymin, ymax);
- gmt_world_map = (fabs (project_info.w - project_info.e) == 360.0);
- gmtdefs.n_lat_nodes = 2;
- gmtdefs.n_lon_nodes = 3; /* > 2 to avoid map-jumps */
- forward = (PFI)merc; inverse = (PFI)imerc;
- outside = (PFI) wesn_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.horizontal = TRUE;
- return (FALSE); /* No need to search for wesn */
- }
-
- int vmerc (cmerid)
- double cmerid; {
- /* Set up a Mercator transformation */
-
- project_info.central_meridian = cmerid;
- return (0);
- }
-
- int merc (lon, lat, x, y)
- double lon, lat, *x, *y; {
- double tmp;
- /* Convert lon/lat to Mercator x/y */
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- *x = (lon - project_info.central_meridian) * D2R * EQ_RAD;
-
- if (fabs (lat) < 90.0) {
- tmp = pow ((1.0 - ECC * sin (lat * D2R)) / (1.0 + ECC * sin (lat * D2R)), 0.5 * ECC);
- *y = d_log (tan ((45.0 + 0.5 * lat) * D2R) * tmp) * EQ_RAD;
- }
- else
- *y = copysign (1.0e100, lat);
- }
-
- int imerc (lon, lat, x, y)
- double *lon, *lat, x, y; {
- int i;
- double t, t_phi, delta, tmp, phi, pw, pow();
-
- /* Convert Mercator x/y to lon/lat */
-
- *lon = (x / EQ_RAD) * R2D + project_info.central_meridian;
-
- t = exp (-y / EQ_RAD);
- t_phi = 90.0 - 2.0 * atan(t) * R2D;
- for (delta = 1.0, i = 0; delta > 0.0000001 && i < 100; i++) {
- tmp = (1.0 - ECC * sin (t_phi * D2R)) / (1.0 + ECC * sin (t_phi * D2R));
- pw = pow (tmp, 0.5 * ECC);
- phi = 90.0 - 2.0 * atan (t * pw) * R2D;
- delta = fabs (fabs (phi) - fabs (t_phi));
- t_phi = phi;
- }
- *lat = phi;
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR CYLIDNDRICAL EQUAL-AREA PROJECTIONS (CYL_EQ)
- */
-
- int map_init_cyleq () {
- double xmin, xmax, ymin, ymax;
-
- gmt_world_map = (fabs (project_info.w - project_info.e) == 360.0);
- vcyleq (project_info.pars[0], project_info.pars[1]);
- cyleq (project_info.w, project_info.s, &xmin, &ymin);
- cyleq (project_info.e, project_info.n, &xmax, &ymax);
- if (project_info.gave_map_width) project_info.pars[2] /= (M_PR_DEG * fabs (project_info.w - project_info.e));
- if (project_info.units_pr_degree) project_info.pars[2] /= M_PR_DEG;
- project_info.x_scale = project_info.y_scale = project_info.pars[2];
- map_setxy (xmin, xmax, ymin, ymax);
- gmtdefs.n_lat_nodes = 2;
- gmtdefs.n_lon_nodes = 3; /* > 2 to avoid map-jumps */
- forward = (PFI)cyleq; inverse = (PFI)icyleq;
- outside = (PFI) wesn_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.horizontal = TRUE;
- return (FALSE); /* No need to search for wesn */
- }
-
- int vcyleq (lon0, slat)
- double lon0, slat; {
- /* Set up a Cylindrical equal-area transformation */
-
- gmt_check_R_J (&lon0);
- project_info.central_meridian = lon0;
- project_info.y_rx = EQ_RAD * D2R * cos (D2R * slat);
- project_info.y_ry = EQ_RAD / cos (D2R * slat);
- return (0);
- }
-
- int cyleq (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to Cylindrical equal-area x/y */
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- *x = (lon - project_info.central_meridian) * project_info.y_rx;
- *y = project_info.y_ry * sind (lat);
- }
-
- int icyleq (lon, lat, x, y)
- double *lon, *lat, x, y; {
-
- /* Convert Cylindrical equal-area x/y to lon/lat */
-
- *lon = (x / project_info.y_rx) + project_info.central_meridian;
- *lat = R2D * d_asin (y / project_info.y_ry);
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR CYLINDRICAL EQIDISTANT PROJECTION (CYL_EQDIST)
- */
-
- int map_init_cyleqdist () {
- double xmin, xmax, ymin, ymax;
-
- gmt_world_map = (fabs (project_info.w - project_info.e) == 360.0);
- vcyleqdist (project_info.pars[0]);
- cyleqdist (project_info.w, project_info.s, &xmin, &ymin);
- cyleqdist (project_info.e, project_info.n, &xmax, &ymax);
- if (project_info.gave_map_width) project_info.pars[1] /= (M_PR_DEG * fabs (project_info.w - project_info.e));
- if (project_info.units_pr_degree) project_info.pars[1] /= M_PR_DEG;
- project_info.x_scale = project_info.y_scale = project_info.pars[1];
- map_setxy (xmin, xmax, ymin, ymax);
- gmtdefs.n_lat_nodes = 2;
- gmtdefs.n_lon_nodes = 3; /* > 2 to avoid map-jumps */
- forward = (PFI)cyleqdist; inverse = (PFI)icyleqdist;
- outside = (PFI) wesn_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.horizontal = TRUE;
- return (FALSE); /* No need to search for wesn */
- }
-
- int vcyleqdist (lon0)
- double lon0; {
- /* Set up a Cylindrical equidistant transformation */
-
- gmt_check_R_J (&lon0);
- project_info.central_meridian = lon0;
- return (0);
- }
-
- int cyleqdist (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to Cylindrical equidistant x/y */
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- *x = (lon - project_info.central_meridian) * D2R * EQ_RAD;
- *y = lat * D2R * EQ_RAD;
- }
-
- int icyleqdist (lon, lat, x, y)
- double *lon, *lat, x, y; {
-
- /* Convert Cylindrical equal-area x/y to lon/lat */
-
- *lon = R2D * (x / EQ_RAD) + project_info.central_meridian;
- *lat = R2D * y / EQ_RAD;
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE POLAR STEREOGRAPHIC PROJECTION
- */
-
- int map_init_stereo () {
- BOOLEAN search;
- double xmin, xmax, ymin, ymax, dummy, radius;
-
- if (fabs (project_info.pars[1]) == 90.0) { /* Polar aspect */
- project_info.polar = TRUE;
- forward = (PFI)plrs; inverse = (PFI)iplrs;
- vplrs (project_info.pars[0], project_info.pars[1]);
- if (project_info.units_pr_degree) {
- plrs (project_info.pars[0], project_info.pars[3], &dummy, &radius);
- project_info.x_scale = project_info.y_scale = fabs (project_info.pars[2] / radius);
- }
- else
- project_info.x_scale = project_info.y_scale = project_info.pars[2];
- }
- else {
- vstereo (project_info.pars[0], project_info.pars[1]);
- project_info.polar = FALSE;
- forward = (project_info.pole == 0.0) ? (PFI)stereo2 : (PFI)stereo1;
- inverse = (PFI)istereo;
- if (project_info.units_pr_degree) {
- vplrs (0.0, 90.0);
- plrs (0.0, fabs(project_info.pars[3]), &dummy, &radius);
- project_info.x_scale = project_info.y_scale = fabs (project_info.pars[2] / radius);
- }
- else
- project_info.x_scale = project_info.y_scale = project_info.pars[2];
-
- vstereo (project_info.pars[0], project_info.pars[1]);
- }
-
- if (!project_info.region) { /* Rectangular box given */
- (*forward) (project_info.w, project_info.s, &xmin, &ymin);
- (*forward) (project_info.e, project_info.n, &xmax, &ymax);
-
- outside = (PFI) rect_outside2;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.check_side = !gmtdefs.oblique_anotation;
- frame_info.horizontal = (fabs (project_info.pars[1]) < 30.0 && fabs (project_info.n - project_info.s) < 30.0);
- search = TRUE;
- }
- else {
- if (fabs (project_info.pars[1]) == 90.0) { /* Polar aspect */
- if (!project_info.north_pole && project_info.s == -90.0) project_info.edge[0] = FALSE;
- if (project_info.north_pole && project_info.n == 90.0) project_info.edge[2] = FALSE;
- if ((fabs (project_info.w - project_info.e) == 360.0 || project_info.w == project_info.e)) project_info.edge[1] = project_info.edge[3] = FALSE;
- outside = (PFI) polar_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- frame_info.horizontal = TRUE;
- gmtdefs.n_lat_nodes = 2;
- xy_search (&xmin, &xmax, &ymin, &ymax);
- }
- else { /* Global view only */
- frame_info.anot_int[0] = frame_info.anot_int[1] = 0.0; /* No annotations for global mode */
- frame_info.frame_int[0] = frame_info.frame_int[1] = 0.0; /* No tickmarks for global mode */
- project_info.w = 0.0;
- project_info.e = 360.0;
- project_info.s = -90.0;
- project_info.n = 90.0;
- xmin = ymin = -2.0 * EQ_RAD;
- xmax = ymax = -xmin;
- outside = (PFI) radial_outside;
- crossing = (PFI) radial_crossing;
- overlap = (PFI) radial_overlap;
- map_clip = (PFI) radial_clip;
- gmtdefs.basemap_type = 1;
- }
- search = FALSE;
- left_edge = (PFD) left_circle;
- right_edge = (PFD) right_circle;
- }
-
- map_setinfo (xmin, xmax, ymin, ymax, project_info.pars[2]);
- project_info.r = 0.5 * project_info.xmax;
-
- return (search);
- }
-
- int vplrs (rlong0, plat)
- double rlong0, plat; {
- /* Set up a polar stereographic transformation (spherical) */
-
- project_info.central_meridian = rlong0;
- project_info.pole = plat;
- project_info.north_pole = (plat > 0.0);
- project_info.s_c = 2.0 * EQ_RAD * gmtdefs.map_scale_factor / d_sqrt (pow (1.0 + ECC, 1.0 + ECC) * pow (1.0 - ECC, 1.0 - ECC));
- return (0);
- }
-
- int plrs (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to x/y using polar projection */
- double t, rho, lambda, phi;
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
-
- lambda = (lon - project_info.central_meridian) * D2R;
- lat *= D2R;
- phi = (project_info.north_pole) ? lat : -lat;
- t = tan (M_PI_4 - 0.5 * phi) / pow ((1.0 - ECC * sin (phi)) / (1.0 + ECC * sin (phi)), 0.5 * ECC);
- rho = project_info.s_c * t;
-
- *x = rho * sin (lambda);
- *y = (project_info.north_pole) ? -rho * cos (lambda) : rho * cos (lambda);
- }
-
- int iplrs (lon, lat, x, y)
- double *lon, *lat, x, y; {
- int i;
- double rho, temp, tphi, t, delta, phi;
-
- /* Convert polar x/y to lon/lat */
-
- if (x == 0.0 && y == 0.0) {
- *lon = project_info.central_meridian;
- *lat = project_info.pole;
- }
- else {
- *lon = project_info.central_meridian + ((project_info.north_pole) ? d_atan2 (x, -y) : d_atan2 (x, y)) * R2D;
- rho = hypot (x, y);
- t = rho / project_info.s_c;
- tphi = M_PI_2 - 2.0 * atan (t);
- delta = 1.0;
- for (i = 0; i < 100 && delta > 1.0e-8; i++) {
- temp = (1.0 - ECC * sin (tphi)) / (1.0 + ECC * sin (tphi));
- phi = M_PI_2 - 2.0 * atan (t * pow (temp, 0.5 * ECC));
- delta = fabs (fabs (tphi) - fabs (phi));
- tphi = phi;
- }
- *lat = (project_info.north_pole) ? phi * R2D : -phi * R2D;
- }
- }
-
- /* STEREO: For general oblique view */
-
- int vstereo (rlong0, plat)
- double rlong0, plat; {
- double s, m1, z1;
- /* Set up a polar stereographic transformation (spherical) */
-
- project_info.central_meridian = rlong0;
- project_info.pole = plat;
- project_info.sinp = sind (plat);
- project_info.cosp = cosd (plat);
- project_info.north_pole = (plat > 0.0);
- m1 = cosd (plat) / sqrt (1.0 - ECC2 * pow (sind (plat), 2.0));
- s = sind (plat);
- z1 = 2.0 * atan (sqrt (pow ((1.0 - ECC * s) / (1.0 + ECC * s), ECC) * (1.0 + s) / (1.0 - s))) - M_PI_2;
- project_info.s_cosz1 = cos (z1);
- project_info.s_sinz1 = sin (z1);
- project_info.s_c = 2.0 * EQ_RAD * gmtdefs.map_scale_factor * m1;
- return (0);
- }
-
- int stereo1 (lon, lat, x, y)
- double lon, lat, *x, *y; {
- double dlon, s, z, A;
-
- /* Convert lon/lat to x/y using stereographic projection, oblique view */
-
- dlon = lon - project_info.central_meridian;
- s = sind (lat);
- z = 2.0 * atan (sqrt (pow ((1.0 - ECC * s) / (1.0 + ECC * s), ECC) * (1.0 + s) / (1.0 - s))) - M_PI_2;
- A = project_info.s_c / (project_info.s_cosz1 * (1.0 + project_info.s_sinz1 * sin (z) +
- project_info.s_cosz1 * cos (z) * cosd (dlon)));
- *x = A * cos (z) * sind (dlon);
- *y = A * (project_info.s_cosz1 * sin (z) - project_info.s_sinz1 * cos (z) * cosd (dlon));
- }
-
- int istereo (lon, lat, x, y)
- double *lon, *lat, x, y; {
- int i;
- double rho, s, z, delta, ce, tphi, phi;
-
- if (x == 0.0 && y == 0.0) {
- *lon = project_info.central_meridian;
- *lat = project_info.pole;
- }
- else {
- rho = hypot (x, y);
- ce = 2.0 * atan (rho * project_info.s_cosz1 / project_info.s_c);
- z = d_asin (cos (ce) * project_info.s_sinz1 + (y * sin (ce) * project_info.s_cosz1 / rho));
- delta = 1.0;
- s = sin (z);
- tphi = 2.0 * atan (tan (M_PI_4 + 0.5 * z) * pow ((1.0 + ECC * s)/(1.0 - ECC * s), 0.5 * ECC)) - M_PI_2;
- for (i = 0; i < 100 && delta > 1.0e-8; i++) {
- s = sin (tphi);
- phi = 2.0 * atan (tan (M_PI_4 + 0.5 * z) * pow ((1.0 + ECC * s)/(1.0 - ECC * s), 0.5 * ECC)) - M_PI_2;
- delta = fabs (fabs (tphi) - fabs (phi));
- tphi = phi;
- }
- *lat = phi * R2D;
- *lon = project_info.central_meridian +
- R2D * atan (x * sin (ce) / (rho * project_info.s_cosz1 * cos (ce) -
- y * project_info.s_sinz1 * sin (ce)));
- }
- }
-
- /* For equatorial view */
-
- int stereo2 (lon, lat, x, y)
- double lon, lat, *x, *y; {
- double dlon, s, z, A;
-
- /* Convert lon/lat to x/y using stereographic projection, equatorial view */
-
- dlon = lon - project_info.central_meridian;
- s = sind (lat);
- z = 2.0 * atan (sqrt (pow ((1.0 - ECC * s) / (1.0 + ECC * s), ECC) * (1.0 + s) / (1.0 - s))) - M_PI_2;
- A = project_info.s_c / (project_info.s_cosz1 * (1.0 + project_info.s_sinz1 * sin (z) +
- project_info.s_cosz1 * cos (z) * cosd (dlon)));
- *x = A * cos (z) * sind (dlon);
- *y = A * sin (z);
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE LAMBERT CONFORMAL CONIC PROJECTION
- */
-
- int map_init_lambert () {
- BOOLEAN search;
- double xmin, xmax, ymin, ymax;
-
- vlamb (project_info.pars[0], project_info.pars[1], project_info.pars[2], project_info.pars[3]);
- if (project_info.units_pr_degree) project_info.pars[4] /= M_PR_DEG;
- project_info.x_scale = project_info.y_scale = project_info.pars[4];
- forward = (PFI)lamb; inverse = (PFI)ilamb;
-
- if (!project_info.region) { /* Rectangular box given*/
- lamb (project_info.w, project_info.s, &xmin, &ymin);
- lamb (project_info.e, project_info.n, &xmax, &ymax);
-
- outside = (PFI) rect_outside;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- search = TRUE;
- frame_info.check_side = !gmtdefs.oblique_anotation;
- }
- else {
- xy_search (&xmin, &xmax, &ymin, &ymax);
- outside = (PFI) wesn_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- left_edge = (PFD) left_lambert;
- right_edge = (PFD) right_lambert;
- search = FALSE;
- }
- map_setinfo (xmin, xmax, ymin, ymax, project_info.pars[4]);
- gmtdefs.n_lat_nodes = 2;
- frame_info.horizontal = TRUE;
- return (search);
- }
-
- int vlamb (rlong0, rlat0, pha, phb)
- double rlong0, rlat0, pha, phb; {
- /* Set up a Lambert Conformal Conic projection */
-
- double t_pha, m_pha, t_phb, m_phb, t_rlat0;
-
- project_info.north_pole = (rlat0 > 0.0);
- project_info.pole = (project_info.north_pole) ? 90.0 : -90.0;
- pha *= D2R;
- phb *= D2R;
-
- t_pha = tan (45.0 * D2R - 0.5 * pha) / pow ((1.0 - ECC *
- sin (pha)) / (1.0 + ECC * sin (pha)), 0.5 * ECC);
- m_pha = cos (pha) / d_sqrt (1.0 - ECC2
- * pow (sin (pha), 2.0));
- t_phb = tan (45.0 * D2R - 0.5 * phb) / pow ((1.0 - ECC *
- sin (phb)) / (1.0 + ECC * sin (phb)), 0.5 * ECC);
- m_phb = cos (phb) / d_sqrt (1.0 - ECC2
- * pow (sin (phb), 2.0));
- t_rlat0 = tan (45.0 * D2R - 0.5 * rlat0 * D2R) /
- pow ((1.0 - ECC * sin (rlat0 * D2R)) /
- (1.0 + ECC * sin (rlat0 * D2R)), 0.5 * ECC);
-
- if(pha != phb)
- project_info.l_N = (d_log(m_pha) - d_log(m_phb))/(d_log(t_pha) - d_log(t_phb));
- else
- project_info.l_N = sin(pha);
-
- project_info.l_F = m_pha/(project_info.l_N * pow(t_pha,project_info.l_N));
- project_info.central_meridian = rlong0;
- project_info.l_rho0 = EQ_RAD * project_info.l_F * pow(t_rlat0,project_info.l_N);
- }
-
- int lamb (lon, lat, x, y)
- double lon, lat, *x, *y; {
- double rho,theta,hold1,hold2,hold3;
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- lat *= D2R;
-
- hold2 = pow (((1.0 - ECC * sin (lat)) / (1.0 + ECC * sin (lat))), 0.5 * ECC);
- hold3 = tan (45.0 * D2R - 0.5 * lat);
- if (fabs (hold3) < SMALL) hold3 = 0.0;
- hold1 = (hold3 == 0.0) ? 0.0 : pow (hold3 / hold2, project_info.l_N);
- rho = EQ_RAD * project_info.l_F * hold1;
- theta = project_info.l_N * (lon - project_info.central_meridian) * D2R;
-
- *x = rho * sin (theta);
- *y = project_info.l_rho0 - rho * cos (theta);
- }
-
- int ilamb (lon, lat, x, y)
- double *lon, *lat, x, y; {
- int i;
- double theta, temp, rho, t, tphi, phi, delta;
-
- theta = atan (x / (project_info.l_rho0 - y));
- *lon = (theta /project_info.l_N) * R2D + project_info.central_meridian;
-
- temp = x * x + (project_info.l_rho0 - y) * (project_info.l_rho0 - y);
- rho = copysign (d_sqrt (temp), project_info.l_N);
- t = pow ((rho / (EQ_RAD * project_info.l_F)), 1./project_info.l_N);
- tphi = 90.0 * D2R - 2.0 * atan (t);
- delta = 1.0;
- for (i = 0; i < 100 && delta > 1.0e-8; i++) {
- temp = (1.0 - ECC * sin (tphi)) / (1.0 + ECC * sin (tphi));
- phi = 90.0 * D2R - 2.0 * atan (t * pow (temp, 0.5 * ECC));
- delta = fabs (fabs (tphi) - fabs (phi));
- tphi = phi;
- }
- *lat = phi * R2D;
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE OBLIQUE MERCATOR PROJECTION
- */
-
- int map_init_oblique () {
- double xmin, xmax, ymin, ymax;
- double o_x, o_y, p_x, p_y, c, az, b_x, b_y, w, e, s, n;
-
- gmt_set_spherical (); /* PW: Force spherical for now */
-
- if (project_info.units_pr_degree) project_info.pars[4] /= M_PR_DEG; /* To get plot-units / m */
-
- o_x = project_info.pars[0]; o_y = project_info.pars[1];
-
- if (project_info.pars[6] == 1.0) { /* Must get correct origin, then get second point */
- p_x = project_info.pars[2]; p_y = project_info.pars[3];
-
- project_info.o_pole_lon = D2R * p_x;
- project_info.o_sin_pole_lat = sind (p_y);
- project_info.o_cos_pole_lat = cosd (p_y);
-
- /* Find azimuth to pole, add 90, and compute second point 10 degrees away */
-
- get_origin (o_x, o_y, p_x, p_y, &o_x, &o_y);
- az = R2D * atan (cosd (p_y) * sind (p_x - o_x) / (cosd (o_y) * sind (p_y) - sind (o_y) * cosd (p_y) * cosd (p_x - o_x))) + 90.0;
- c = 10.0; /* compute point 10 degrees from origin along azimuth */
- b_x = o_x + R2D * atan (sind (c) * sind (az) / (cosd (o_y) * cosd (c) - sind (o_y) * sind (c) * cosd (az)));
- b_y = R2D * d_asin (sind (o_y) * cosd (c) + cosd (o_y) * sind (c) * cosd (az));
- project_info.pars[0] = o_x; project_info.pars[1] = o_y;
- project_info.pars[2] = b_x; project_info.pars[3] = b_y;
- }
- else { /* Just find pole */
- b_x = project_info.pars[2]; b_y = project_info.pars[3];
- get_rotate_pole (o_x, o_y, b_x, b_y);
- }
-
- vmerc (0.0); /* Because origin will have 0 degree oblique */
-
- if (project_info.region) {
- /* Convert oblique wesn in degrees to meters using regular Mercator */
- if (fabs (project_info.w - project_info.e) == 360.0) {
- project_info.w = -180.0;
- project_info.e = 180.0;
- }
- merc (project_info.w, project_info.s, &xmin, &ymin);
- merc (project_info.e, project_info.n, &xmax, &ymax);
- vmerc (0.0);
- project_info.region = FALSE; /* Since wesn was oblique, not geographical wesn */
- }
- else {
- /* wesn is lower left and upper right corners in normal lon/lats */
-
- oblmrc (project_info.w, project_info.s, &xmin, &ymin);
- oblmrc (project_info.e, project_info.n, &xmax, &ymax);
- }
-
- imerc (&w, &s, xmin, ymin); /* Get oblique wesn boundaries */
- imerc (&e, &n, xmax, ymax);
- project_info.x_scale = project_info.y_scale = project_info.pars[4];
- map_setinfo (xmin, xmax, ymin, ymax, project_info.pars[4]);
-
- forward = (PFI)oblmrc; inverse = (PFI)ioblmrc;
- outside = (PFI) rect_outside;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
-
- gmt_world_map = (fabs (fabs (w - e) - 360.0) < SMALL);
- gmtdefs.basemap_type = 1;
- frame_info.check_side = !gmtdefs.oblique_anotation;
-
- return (TRUE);
- }
-
- int oblmrc (lon, lat, x, y)
- double lon, lat;
- double *x, *y; {
- /* Convert a longitude/latitude point to Oblique Mercator coordinates
- * by way of rotation coordinates and then using regular Mercator */
- double tlon, tlat;
-
- pole_rotate_forward (lon, lat, &tlon, &tlat);
- merc (tlon, tlat, x, y);
- }
-
- int ioblmrc (lon, lat, x, y)
- double *lon, *lat;
- double x, y; {
- /* Convert a longitude/latitude point from Oblique Mercator coordinates
- * by way of regular Mercator and then rotate coordinates */
- double tlon, tlat;
-
- imerc (&tlon, &tlat, x, y);
- pole_rotate_inverse (lon, lat, tlon, tlat);
- }
-
- int pole_rotate_forward (lon, lat, tlon, tlat)
- double lon, lat, *tlon, *tlat; {
- /* Given the pole position in project_info, geographical coordinates
- * are computed from oblique coordinates assuming a spherical earth.
- */
-
- double dlon, x, y, test;
-
- lon *= D2R; lat *= D2R;
- dlon = lon - project_info.o_pole_lon;
- test = project_info.o_sin_pole_lat * sin (lat) + project_info.o_cos_pole_lat * cos (lat) * cos (dlon);
- y = d_asin (test);
- x = project_info.o_beta + d_atan2 (cos (lat) * sin (dlon), project_info.o_sin_pole_lat * cos (lat) * cos (dlon)
- - project_info.o_cos_pole_lat * sin (lat));
-
- *tlon = R2D * x; *tlat = R2D * y;
- }
-
- int pole_rotate_inverse (lon, lat, tlon, tlat)
- double *lon, *lat, tlon, tlat; {
- /* Given the pole position in project_info, geographical coordinates
- * are computed from oblique coordinates assuming a spherical earth.
- */
-
- double dlon, x, y, test;
-
- tlon *= D2R; tlat *= D2R;
- dlon = tlon - project_info.o_beta;
- test = project_info.o_sin_pole_lat * sin (tlat) - project_info.o_cos_pole_lat * cos (tlat) * cos (dlon);
- y = d_asin (test);
- x = project_info.o_pole_lon + d_atan2 (cos (tlat) * sin (dlon), project_info.o_sin_pole_lat * cos (tlat) * cos (dlon)
- + project_info.o_cos_pole_lat * sin (tlat));
-
- *lon = R2D * x; *lat = R2D * y;
- }
-
- int get_rotate_pole (lon1, lat1, lon2, lat2)
- double lon1, lat1, lon2, lat2; {
- double plon, plat, beta, dummy;
-
- #ifdef aix
- double aix_cpp_sucks1, aix_cpp_sucks2;
-
- aix_cpp_sucks1 = cosd (lat1) * sind (lat2) * cosd (lon1) - sind (lat1) * cosd (lat2) * cosd (lon2);
- aix_cpp_sucks2 = sind (lat1) * cosd (lat2) * sind (lon2) - cosd (lat1) * sind (lat2) * sind (lon1);
- plon = R2D * d_atan2 (aix_cpp_sucks1, aix_cpp_sucks2);
- #else
- plon = R2D * d_atan2 (cosd (lat1) * sind (lat2) * cosd (lon1) - sind (lat1) * cosd (lat2) * cosd (lon2),
- sind (lat1) * cosd (lat2) * sind (lon2) - cosd (lat1) * sind (lat2) * sind (lon1));
- #endif
-
- plat = R2D * atan (-cosd (plon - lon1) / tand (lat1));
- if (plat < 0.0) {
- plat = -plat;
- plon += 180.0;
- if (plon >= 360.0) plon -= 360.0;
- }
- project_info.o_pole_lon = D2R * plon;
- project_info.o_sin_pole_lat = sind (plat);
- project_info.o_cos_pole_lat = cosd (plat);
- pole_rotate_forward (lon1, lat1, &beta, &dummy);
- project_info.o_beta = -beta * D2R;
- }
-
- int get_origin (lon1, lat1, lon_p, lat_p, lon2, lat2)
- double lon1, lat1, lon_p, lat_p, *lon2, *lat2; {
- double beta, dummy, d, az, c;
-
-
- /* Now find origin that is 90 degrees from pole, let oblique lon=0 go through lon1/lat1 */
- #ifdef aix
- c = cosd (lat_p) * cosd (lat1) * cosd (lon1-lon_p);
- c = R2D * d_acos (sind (lat_p) * sind (lat1) + c);
- #else
- c = R2D * d_acos (sind (lat_p) * sind (lat1) + cosd (lat_p) * cosd (lat1) * cosd (lon1-lon_p));
- #endif
-
- if (c != 90.0) { /* Not true origin */
- d = fabs (90.0 - c);
- az = R2D * d_asin (sind (lon_p-lon1) * cosd (lat_p) / sind (c));
- if (c < 90.0) az += 180.0;
- *lat2 = R2D * d_asin (sind (lat1) * cosd (d) + cosd (lat1) * sind (d) * cosd (az));
- *lon2 = lon1 + R2D * d_atan2 (sind (d) * sind (az), cosd (lat1) * cosd (d) - sind (lat1) * sind (d) * cosd (az));
- if (gmtdefs.verbose) fprintf (stderr, "%s: GMT Warning: Correct projection origin = %lg/%lg\n", gmt_program, *lon2, *lat2);
- }
- else {
- *lon2 = lon1;
- *lat2 = lat1;
- }
-
- pole_rotate_forward (*lon2, *lat2, &beta, &dummy);
-
-
- project_info.o_beta = -beta * D2R;
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE TRANSVERSE MERCATOR PROJECTION (TM)
- */
-
- int map_init_tm () {
- BOOLEAN search;
- double xmin, xmax, ymin, ymax;
-
- vtm (project_info.pars[0]);
- if (project_info.units_pr_degree) project_info.pars[1] /= M_PR_DEG;
- project_info.x_scale = project_info.y_scale = project_info.pars[1];
- forward = (PFI)tm; inverse = (PFI)itm;
-
- if (project_info.region) {
- xy_search (&xmin, &xmax, &ymin, &ymax);
- outside = (PFI) wesn_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- search = FALSE;
- }
- else { /* Find min values */
- (*forward) (project_info.w, project_info.s, &xmin, &ymin);
- (*forward) (project_info.e, project_info.n, &xmax, &ymax);
- outside = (PFI) rect_outside;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.check_side = TRUE;
- search = TRUE;
- }
-
- frame_info.horizontal = TRUE;
- map_setinfo (xmin, xmax, ymin, ymax, project_info.pars[1]);
-
- gmtdefs.basemap_type = 1;
-
- return (search);
- }
-
- int vtm (lon0)
- double lon0; {
- /* Set up an TM projection */
- double e1;
-
- e1 = (1.0 - d_sqrt (1.0 - ECC2)) / (1.0 + d_sqrt (1.0 - ECC2));
- project_info.t_e2 = ECC2 / (1.0 - ECC2);
- project_info.t_c1 = (1.0 - 0.25 * ECC2 - 3.0 * ECC4 / 64.0 - 5.0 * ECC6 / 256.0);
- project_info.t_c2 = (3.0 * ECC2 / 8.0 + 3.0 * ECC4 / 32.0 + 45.0 * ECC6 / 1024.0);
- project_info.t_c3 = (15.0 * ECC4 / 256.0 + 45.0 * ECC6 / 1024.0);
- project_info.t_c4 = (35.0 * ECC6 / 3072.0);
- project_info.t_ic1 = (1.5 * e1 - 27.0 * pow (e1, 3.0) / 32.0);
- project_info.t_ic2 = (21.0 * e1 * e1 / 16.0 - 55.0 * pow (e1, 4.0) / 32.0);
- project_info.t_ic3 = (151.0 * pow (e1, 3.0) / 96.0);
- project_info.t_ic4 = (1097.0 * pow (e1, 4.0) / 512.0);
- project_info.central_meridian = lon0;
- }
-
- int tm (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to TM x/y */
- double N, T, T2, C, A, M, dlon, tan_lat, cos_lat, A2, A3, A5;
-
- dlon = lon - project_info.central_meridian;
- if (fabs (dlon) > 360.0) dlon += copysign (360.0, -dlon);
- if (fabs (dlon) > 180.0) dlon = copysign (360.0 - fabs (dlon), -dlon);
- lat *= D2R;
- M = EQ_RAD * (project_info.t_c1 * lat - project_info.t_c2 * sin (2.0 * lat)
- + project_info.t_c3 * sin (4.0 * lat) - project_info.t_c4 * sin (6.0 * lat));
- if (fabs (lat) == M_PI_2) {
- *x = 0.0;
- *y = gmtdefs.map_scale_factor * M;
- }
- else {
- N = EQ_RAD / d_sqrt (1.0 - ECC2 * pow (sin (lat), 2.0));
- tan_lat = tan (lat);
- cos_lat = cos (lat);
- T = tan_lat * tan_lat;
- T2 = T * T;
- C = project_info.t_e2 * cos_lat * cos_lat;
- A = dlon * D2R * cos_lat;
- A2 = A * A; A3 = A2 * A; A5 = A3 * A2;
- *x = gmtdefs.map_scale_factor * N * (A + (1.0 - T + C) * (A3 * 0.16666666666666666667)
- + (5.0 - 18.0 * T + T2 + 72.0 * C - 58.0 * project_info.t_e2) * (A5 * 0.00833333333333333333));
- A3 *= A; A5 *= A;
- *y = gmtdefs.map_scale_factor * (M + N * tan (lat) * (0.5 * A2 + (5.0 - T + 9.0 * C + 4.0 * C * C) * (A3 * 0.04166666666666666667)
- + (61.0 - 58.0 * T + T2 + 600.0 * C - 330.0 * project_info.t_e2) * (A5 * 0.00138888888888888889)));
- }
- }
-
- int itm (lon, lat, x, y)
- double *lon, *lat, x, y; {
- /* Convert TM x/y to lon/lat */
- double M, mu, phi1, C1, C12, T1, T12, tmp, tmp2, N1, R1, D, D2, D3, D5, cos_phi1, tan_phi1;
-
- M = y / gmtdefs.map_scale_factor;
- mu = M / (EQ_RAD * project_info.t_c1);
- phi1 = mu + project_info.t_ic1 * sin (2.0 * mu) + project_info.t_ic2 * sin (4.0 * mu)
- + project_info.t_ic3 * sin (6.0 * mu) + project_info.t_ic4 * sin (8.0 * mu);
- cos_phi1 = cos (phi1);
- tan_phi1 = tan (phi1);
- C1 = project_info.t_e2 * cos_phi1 * cos_phi1;
- C12 = C1 * C1;
- T1 = tan_phi1 * tan_phi1;
- T12 = T1 * T1;
- tmp = 1.0 - ECC2 * (1.0 - cos_phi1 * cos_phi1);
- tmp2 = d_sqrt (tmp);
- N1 = EQ_RAD / tmp2;
- R1 = EQ_RAD * (1.0 - ECC2) / (tmp * tmp2);
- D = x / (N1 * gmtdefs.map_scale_factor);
- D2 = D * D; D3 = D2 * D; D5 = D3 * D2;
-
- *lon = project_info.central_meridian + R2D * (D - (1.0 + 2.0 * T1 + C1) * (D3 * 0.16666666666666666667)
- + (5.0 - 2.0 * C1 + 28.0 * T1 - 3.0 * C12 + 8.0 * project_info.t_e2 + 24.0 * T12)
- * (D5 * 0.00833333333333333333)) / cos (phi1);
- D3 *= D; D5 *= D;
- *lat = phi1 - (N1 * tan (phi1) / R1) * (0.5 * D2 -
- (5.0 + 3.0 * T1 + 10.0 * C1 - 4.0 * C12 - 9.0 * project_info.t_e2) * (D3 * 0.04166666666666666667)
- + (61.0 + 90.0 * T1 + 298 * C1 + 45.0 * T12 - 252.0 * project_info.t_e2 - 3.0 * C12) * (D5 * 0.00138888888888888889));
- (*lat) *= R2D;
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE UNIVERSAL TRANSVERSE MERCATOR PROJECTION (UTM)
- */
-
- int map_init_utm () {
- BOOLEAN search;
- double xmin, xmax, ymin, ymax, lon0;
-
- lon0 = 180.0 + 6.0 * project_info.pars[0] - 3.0;
- if (lon0 >= 360.0) lon0 -= 360.0;
- vtm (lon0); /* Central meridian for this zone */
- if (project_info.units_pr_degree) project_info.pars[1] /= M_PR_DEG;
- project_info.x_scale = project_info.y_scale = project_info.pars[1];
- forward = (PFI)utm; inverse = (PFI)iutm;
-
- if (fabs (project_info.w - project_info.e) > 360.0) { /* -R in UTM meters */
- iutm (&project_info.w, &project_info.s, project_info.w, project_info.s);
- iutm (&project_info.e, &project_info.n, project_info.e, project_info.n);
- project_info.region = FALSE;
- }
- if (project_info.region) {
- xy_search (&xmin, &xmax, &ymin, &ymax);
- outside = (PFI) wesn_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- search = FALSE;
- }
- else {
- utm (project_info.w, project_info.s, &xmin, &ymin);
- utm (project_info.e, project_info.n, &xmax, &ymax);
- outside = (PFI) rect_outside;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.check_side = TRUE;
- search = TRUE;
- }
-
- frame_info.horizontal = TRUE;
- map_setxy (xmin, xmax, ymin, ymax);
-
- gmtdefs.basemap_type = 1;
-
- return (search);
- }
-
- int utm (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to UTM x/y */
-
- if (lon < 0.0) lon += 360.0;
- tm (lon, lat, x, y);
- (*x) += 500000.0;
- if (!project_info.north_pole) (*y) += 10000000.0; /* For S hemisphere, add 10^6 m */
- }
-
- int iutm (lon, lat, x, y)
- double *lon, *lat, x, y; {
- /* Convert UTM x/y to lon/lat */
-
- x -= 500000.0;
- if (!project_info.north_pole) y -= 10000000.0;
- itm (lon, lat, x, y);
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE LAMBERT AZIMUTHAL EQUAL AREA PROJECTION
- */
-
- int map_init_lambeq () {
- BOOLEAN search;
- double xmin, xmax, ymin, ymax, dummy, radius;
-
- gmt_set_spherical (); /* PW: Force spherical for now */
-
- forward = (PFI)lambeq; inverse = (PFI)ilambeq;
- if (project_info.units_pr_degree) {
- vlambeq (0.0, 90.0);
- lambeq (0.0, fabs(project_info.pars[3]), &dummy, &radius);
- project_info.x_scale = project_info.y_scale = fabs (project_info.pars[2] / radius);
- }
- else
- project_info.x_scale = project_info.y_scale = project_info.pars[2];
- vlambeq (project_info.pars[0], project_info.pars[1]);
-
- if (fabs (project_info.pars[1]) == 90.0) {
- project_info.polar = TRUE;
- project_info.north_pole = (project_info.pars[1] == 90.0);
- }
- if (!project_info.region) { /* Rectangular box given */
- (*forward) (project_info.w, project_info.s, &xmin, &ymin);
- (*forward) (project_info.e, project_info.n, &xmax, &ymax);
-
- outside = (PFI) rect_outside2;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.check_side = !gmtdefs.oblique_anotation;
- frame_info.horizontal = (fabs (project_info.pars[1]) < 30.0 && fabs (project_info.n - project_info.s) < 30.0);
- search = TRUE;
- }
- else {
- if (project_info.polar && (project_info.n - project_info.s) < 180.0) { /* Polar aspect */
- if (!project_info.north_pole && project_info.s == -90.0) project_info.edge[0] = FALSE;
- if (project_info.north_pole && project_info.n == 90.0) project_info.edge[2] = FALSE;
- if ((fabs (project_info.w - project_info.e) == 360.0 || project_info.w == project_info.e)) project_info.edge[1] = project_info.edge[3] = FALSE;
- outside = (PFI) polar_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- frame_info.horizontal = TRUE;
- gmtdefs.n_lat_nodes = 2;
- xy_search (&xmin, &xmax, &ymin, &ymax);
- }
- else { /* Global view only */
- frame_info.anot_int[0] = frame_info.anot_int[1] = 0.0; /* No annotations for global mode */
- frame_info.frame_int[0] = frame_info.frame_int[1] = 0.0; /* No tickmarks for global mode */
- project_info.w = 0.0;
- project_info.e = 360.0;
- project_info.s = -90.0;
- project_info.n = 90.0;
- xmin = ymin = -d_sqrt (2.0) * EQ_RAD;
- xmax = ymax = -xmin;
- outside = (PFI) radial_outside;
- crossing = (PFI) radial_crossing;
- overlap = (PFI) radial_overlap;
- map_clip = (PFI) radial_clip;
- gmtdefs.basemap_type = 1;
- }
- search = FALSE;
- left_edge = (PFD) left_circle;
- right_edge = (PFD) right_circle;
- }
-
- map_setinfo (xmin, xmax, ymin, ymax, project_info.pars[2]);
- project_info.r = 0.5 * project_info.xmax;
-
- return (search);
- }
-
- int vlambeq (lon0, lat0)
- double lon0, lat0; {
- /* Set up Lambert Azimuthal Equal-Area projection */
-
- project_info.central_meridian = lon0;
- project_info.pole = lat0;
- project_info.sinp = sin (lat0 * D2R);
- project_info.cosp = cos (lat0 * D2R);
- }
-
- int lambeq (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to Lambert Azimuthal Equal-Area x/y */
- double k, dlon, tmp;
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- dlon = (lon - project_info.central_meridian) * D2R;
- lat *= D2R;
-
- tmp = 1.0 + project_info.sinp * sin (lat) + project_info.cosp * cos (lat) * cos (dlon);
- if (tmp > 0.0) {
- k = d_sqrt (2.0 / tmp);
- *x = EQ_RAD * k * cos (lat) * sin (dlon);
- *y = EQ_RAD * k * (project_info.cosp * sin (lat) - project_info.sinp * cos (lat) * cos (dlon));
- }
- else
- *x = *y = -MAXDOUBLE;
- }
-
- int ilambeq (lon, lat, x, y)
- double *lon, *lat, x, y; {
- /* Convert Lambert Azimuthal Equal-Area x/yto lon/lat */
- double rho, c;
-
- rho = hypot (x, y);
- c = 2.0 * d_asin (0.5 * rho / EQ_RAD);
-
- if (rho == 0.0) {
- *lat = project_info.pole;
- *lon = project_info.central_meridian;
- }
- else {
- *lat = d_asin (cos (c) * project_info.sinp + (y * sin (c) * project_info.cosp / rho)) * R2D;
- if (project_info.pole == 90.0)
- *lon = project_info.central_meridian + R2D * d_atan2 (x, -y);
- else if (project_info.pole == -90.0)
- *lon = project_info.central_meridian + R2D * d_atan2 (x, y);
- else
- *lon = project_info.central_meridian +
- R2D * d_atan2 (x * sin (c), (rho * project_info.cosp * cos (c) - y * project_info.sinp * sin (c)));
- }
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE ORTHOGRAPHIC PROJECTION
- */
-
- int map_init_ortho () {
- BOOLEAN search;
- double xmin, xmax, ymin, ymax, dummy, radius;
-
- gmt_set_spherical (); /* PW: Force spherical for now */
-
- if (project_info.units_pr_degree) {
- vortho (0.0, 90.0);
- ortho (0.0, fabs(project_info.pars[3]), &dummy, &radius);
- project_info.x_scale = project_info.y_scale = fabs (project_info.pars[2] / radius);
- }
- else
- project_info.x_scale = project_info.y_scale = project_info.pars[2];
-
- vortho (project_info.pars[0], project_info.pars[1]);
- forward = (PFI)ortho; inverse = (PFI)iortho;
-
- if (fabs (project_info.pars[1]) == 90.0) {
- project_info.polar = TRUE;
- project_info.north_pole = (project_info.pars[1] == 90.0);
- }
-
- if (!project_info.region) { /* Rectangular box given */
- (*forward) (project_info.w, project_info.s, &xmin, &ymin);
- (*forward) (project_info.e, project_info.n, &xmax, &ymax);
-
- outside = (PFI) rect_outside2;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.check_side = !gmtdefs.oblique_anotation;
- frame_info.horizontal = (fabs (project_info.pars[1]) < 30.0 && fabs (project_info.n - project_info.s) < 30.0);
- search = TRUE;
- }
- else {
- if (project_info.polar) { /* Polar aspect */
- if (!project_info.north_pole && project_info.s == -90.0) project_info.edge[0] = FALSE;
- if (project_info.north_pole && project_info.n == 90.0) project_info.edge[2] = FALSE;
- if ((fabs (project_info.w - project_info.e) == 360.0 || project_info.w == project_info.e)) project_info.edge[1] = project_info.edge[3] = FALSE;
- outside = (PFI) polar_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- frame_info.horizontal = TRUE;
- gmtdefs.n_lat_nodes = 2;
- xy_search (&xmin, &xmax, &ymin, &ymax);
- }
- else { /* Global view only */
- frame_info.anot_int[0] = frame_info.anot_int[1] = 0.0; /* No annotations for global mode */
- frame_info.frame_int[0] = frame_info.frame_int[1] = 0.0; /* No tickmarks for global mode */
- project_info.w = 0.0;
- project_info.e = 360.0;
- project_info.s = -90.0;
- project_info.n = 90.0;
- xmin = ymin = -EQ_RAD;
- xmax = ymax = -xmin;
- outside = (PFI) radial_outside;
- crossing = (PFI) radial_crossing;
- overlap = (PFI) radial_overlap;
- map_clip = (PFI) radial_clip;
- gmtdefs.basemap_type = 1;
- }
- search = FALSE;
- left_edge = (PFD) left_circle;
- right_edge = (PFD) right_circle;
- }
-
- map_setinfo (xmin, xmax, ymin, ymax, project_info.pars[2]);
- project_info.r = 0.5 * project_info.xmax;
-
- return (search);
- }
-
- int vortho (lon0, lat0)
- double lon0, lat0; {
- /* Set up Orthographic projection */
-
- project_info.central_meridian = lon0;
- project_info.pole = lat0;
- project_info.sinp = sin (lat0 * D2R);
- project_info.cosp = cos (lat0 * D2R);
- }
-
- int ortho (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to Orthographic x/y */
- double dlon;
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- dlon = (lon - project_info.central_meridian) * D2R;
- lat *= D2R;
-
- *x = EQ_RAD * cos (lat) * sin (dlon);
- *y = EQ_RAD * (project_info.cosp * sin (lat) - project_info.sinp * cos (lat) * cos (dlon));
- }
-
- int iortho (lon, lat, x, y)
- double *lon, *lat, x, y; {
- /* Convert Orthographic x/yto lon/lat */
- double rho, c;
-
- rho = hypot (x, y);
- c = d_asin (rho / EQ_RAD);
-
- if (rho == 0.0) {
- *lat = project_info.pole;
- *lon = project_info.central_meridian;
- }
- else {
- *lat = d_asin (cos (c) * project_info.sinp + (y * sin (c) * project_info.cosp / rho)) * R2D;
- if (project_info.pole == 90.0)
- *lon = project_info.central_meridian + R2D * d_atan2 (x, -y);
- else if (project_info.pole == -90.0)
- *lon = project_info.central_meridian + R2D * d_atan2 (x, y);
- else
- *lon = project_info.central_meridian +
- R2D * d_atan2 (x * sin (c), (rho * project_info.cosp * cos (c) - y * project_info.sinp * sin (c)));
- }
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE AZIMUTHAL EQUIDISTANT PROJECTION
- */
-
- int map_init_azeqdist () {
- BOOLEAN search;
- double xmin, xmax, ymin, ymax, dummy, radius;
-
- gmt_set_spherical (); /* PW: Force spherical for now */
-
- if (project_info.units_pr_degree) {
- vazeqdist (0.0, 90.0);
- azeqdist (0.0, project_info.pars[3], &dummy, &radius);
- if (radius == 0.0) radius = M_PI * EQ_RAD;
- project_info.x_scale = project_info.y_scale = fabs (project_info.pars[2] / radius);
- }
- else
- project_info.x_scale = project_info.y_scale = project_info.pars[2];
-
- vazeqdist (project_info.pars[0], project_info.pars[1]);
- forward = (PFI)azeqdist; inverse = (PFI)iazeqdist;
-
- if (fabs (project_info.pars[1]) == 90.0) {
- project_info.polar = TRUE;
- project_info.north_pole = (project_info.pars[1] == 90.0);
- }
-
- if (!project_info.region) { /* Rectangular box given */
- (*forward) (project_info.w, project_info.s, &xmin, &ymin);
- (*forward) (project_info.e, project_info.n, &xmax, &ymax);
-
- outside = (PFI) rect_outside;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.check_side = !gmtdefs.oblique_anotation;
- frame_info.horizontal = (fabs (project_info.pars[1]) < 60.0 && fabs (project_info.n - project_info.s) < 30.0);
- search = TRUE;
- }
- else {
- if (project_info.polar && (project_info.n - project_info.s) < 180.0) { /* Polar aspect */
- if (!project_info.north_pole && project_info.s == -90.0) project_info.edge[0] = FALSE;
- if (project_info.north_pole && project_info.n == 90.0) project_info.edge[2] = FALSE;
- if ((fabs (project_info.w - project_info.e) == 360.0 || project_info.w == project_info.e)) project_info.edge[1] = project_info.edge[3] = FALSE;
- outside = (PFI) polar_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- frame_info.horizontal = TRUE;
- gmtdefs.n_lat_nodes = 2;
- xy_search (&xmin, &xmax, &ymin, &ymax);
- }
- else { /* Global view only, force wesn = 0/360/-90/90 */
- frame_info.anot_int[0] = frame_info.anot_int[1] = 0.0; /* No annotations for global mode */
- frame_info.frame_int[0] = frame_info.frame_int[1] = 0.0; /* No tickmarks for global mode */
- project_info.w = 0.0;
- project_info.e = 360.0;
- project_info.s = -90.0;
- project_info.n = 90.0;
- xmin = ymin = -M_PI * EQ_RAD;
- xmax = ymax = -xmin;
- outside = (PFI) eqdist_outside;
- crossing = (PFI) eqdist_crossing;
- overlap = (PFI) radial_overlap;
- map_clip = (PFI) radial_clip;
- gmtdefs.basemap_type = 1;
- }
- search = FALSE;
- left_edge = (PFD) left_circle;
- right_edge = (PFD) right_circle;
- }
-
- map_setinfo (xmin, xmax, ymin, ymax, project_info.pars[2]);
- project_info.r = 0.5 * project_info.xmax;
-
- return (search);
- }
-
- int vazeqdist (lon0, lat0)
- double lon0, lat0; {
- /* Set up azimuthal equidistant projection */
-
- project_info.central_meridian = lon0;
- project_info.pole = lat0;
- project_info.sinp = sin (lat0 * D2R);
- project_info.cosp = cos (lat0 * D2R);
- }
-
- int azeqdist (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to azimuthal equidistant x/y */
- double k, dlon, cc, c, clat, clon, slat;
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- dlon = (lon - project_info.central_meridian) * D2R;
- lat *= D2R;
- slat = sin (lat); clat = cos (lat); clon = cos (dlon);
-
- cc = project_info.sinp * slat + project_info.cosp * clat * clon;
- if (fabs (cc) >= 1.0)
- *x = *y = 0.0;
- else {
- c = d_acos (cc);
- k = EQ_RAD * c / sin (c);
- *x = k * clat * sin (dlon);
- *y = k * (project_info.cosp * slat - project_info.sinp * clat * clon);
- }
- }
-
- int iazeqdist (lon, lat, x, y)
- double *lon, *lat, x, y; {
- /* Convert azimuthal equidistant x/yto lon/lat */
- double rho, c, sin_c, cos_c;
-
- rho = hypot (x, y);
-
- if (rho == 0.0) {
- *lat = project_info.pole;
- *lon = project_info.central_meridian;
- }
- else {
- c = rho / EQ_RAD;
- sin_c = sin (c);
- cos_c = cos (c);
- *lat = d_asin (cos_c * project_info.sinp + (y * sin_c * project_info.cosp / rho)) * R2D;
- if (project_info.pole == 90.0)
- *lon = project_info.central_meridian + R2D * d_atan2 (x, -y);
- else if (project_info.pole == -90.0)
- *lon = project_info.central_meridian + R2D * d_atan2 (x, y);
- else
- *lon = project_info.central_meridian +
- R2D * d_atan2 (x * sin_c, (rho * project_info.cosp * cos_c - y * project_info.sinp * sin_c));
- if ((*lon) <= -180) (*lon) += 360.0;
- }
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE MOLLWEIDE EQUAL AREA PROJECTION
- */
-
- int map_init_mollweide () {
- int search;
- double xmin, xmax, ymin, ymax, y, dummy;
-
- gmt_set_spherical (); /* PW: Force spherical for now */
-
- if (project_info.pars[0] < 0.0) project_info.pars[0] += 360.0;
- gmt_world_map = (fabs (project_info.w - project_info.e) == 360.0);
- if (project_info.gave_map_width) project_info.pars[1] /= (M_PR_DEG * fabs (project_info.w - project_info.e));
- if (project_info.units_pr_degree) project_info.pars[1] /= M_PR_DEG;
- project_info.x_scale = project_info.y_scale = M_PI * project_info.pars[1] / sqrt (8.0);
- vmollweide (project_info.pars[0], project_info.pars[1]);
- if (project_info.s == -90.0) project_info.edge[0] = FALSE;
- if (project_info.n == 90.0) project_info.edge[2] = FALSE;
-
- if (project_info.region) {
- y = (project_info.s * project_info.n <= 0.0) ? 0.0 : MIN (fabs(project_info.s), fabs(project_info.n));
- mollweide (project_info.w, y, &xmin, &dummy);
- mollweide (project_info.e, y, &xmax, &dummy);
- mollweide (project_info.central_meridian, project_info.s, &dummy, &ymin);
- mollweide (project_info.central_meridian, project_info.n, &dummy, &ymax);
- outside = (PFI) wesn_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- left_edge = (PFD) left_ellipse;
- right_edge = (PFD) right_ellipse;
- frame_info.horizontal = 2;
- project_info.polar = TRUE;
- search = FALSE;
- }
- else {
- mollweide (project_info.w, project_info.s, &xmin, &ymin);
- mollweide (project_info.e, project_info.n, &xmax, &ymax);
- outside = (PFI) rect_outside;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.check_side = TRUE;
- search = TRUE;
- }
- map_setxy (xmin, xmax, ymin, ymax);
- forward = (PFI)mollweide; inverse = (PFI)imollweide;
- gmtdefs.basemap_type = 1;
- return (search);
- }
-
- int vmollweide (lon0, scale)
- double lon0, scale; {
- /* Set up Mollweide Equal-Area projection */
-
- gmt_check_R_J (&lon0);
- project_info.central_meridian = lon0;
- project_info.w_x = EQ_RAD * d_sqrt (8.0) / M_PI;
- project_info.w_y = EQ_RAD * M_SQRT2;
- project_info.w_r = 0.25 * (scale * M_PR_DEG * 360.0); /* = Half the minor axis */
- }
-
- int mollweide (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to Mollweide Equal-Area x/y */
- int i = 0;
- double phi, delta;
-
- if (fabs(lat) == 90.0) { /* Special case */
- *x = 0.0;
- *y = copysign (project_info.w_y, lat);
- return;
- }
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- lat *= D2R;
-
- phi = lat;
- do {
- i++;
- delta = -(phi + sin (phi) - M_PI * sin (lat)) / (1.0 + cos (phi));
- phi += delta;
- }
- while (fabs (delta) > 1.0e-7 && i < 100);
- phi *= 0.5;
- *x = project_info.w_x * (lon - project_info.central_meridian) * D2R * cos (phi);
- *y = project_info.w_y * sin (phi);
- }
-
- int imollweide (lon, lat, x, y)
- double *lon, *lat, x, y; {
- /* Convert Mollweide Equal-Area x/y to lon/lat */
- double phi, phi2;
-
- phi = d_asin (y / project_info.w_y);
- phi2 = 2.0 * phi;
- *lat = R2D * d_asin ((phi2 + sin (phi2)) / M_PI);
- *lon = project_info.central_meridian + R2D * (x / (project_info.w_x * cos (phi)));
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE HAMMER-AITOFF EQUAL AREA PROJECTION
- */
-
- int map_init_hammer () {
- int search;
- double xmin, xmax, ymin, ymax, x, y, dummy;
-
- gmt_set_spherical (); /* PW: Force spherical for now */
-
- if (project_info.pars[0] < 0.0) project_info.pars[0] += 360.0;
- gmt_world_map = (fabs (project_info.w - project_info.e) == 360.0);
- if (project_info.gave_map_width) project_info.pars[1] /= (M_PR_DEG * fabs (project_info.w - project_info.e));
- if (project_info.units_pr_degree) project_info.pars[1] /= M_PR_DEG;
- project_info.x_scale = project_info.y_scale = 0.5 * M_PI * project_info.pars[1] / M_SQRT2;
- vhammer (project_info.pars[0], project_info.pars[1]);
- if (project_info.s == -90.0) project_info.edge[0] = FALSE;
- if (project_info.n == 90.0) project_info.edge[2] = FALSE;
-
- if (project_info.region) {
- y = (project_info.s * project_info.n <= 0.0) ? 0.0 : MIN (fabs(project_info.s), fabs(project_info.n));
- x = (fabs (project_info.w - project_info.central_meridian) > fabs (project_info.e - project_info.central_meridian)) ? project_info.w : project_info.e;
- hammer (project_info.w, y, &xmin, &dummy);
- hammer (project_info.e, y, &xmax, &dummy);
- hammer (x, project_info.s, &dummy, &ymin);
- hammer (x, project_info.n, &dummy, &ymax);
- outside = (PFI) wesn_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- left_edge = (PFD) left_ellipse;
- right_edge = (PFD) right_ellipse;
- frame_info.horizontal = 2;
- project_info.polar = TRUE;
- search = FALSE;
- }
- else {
- hammer (project_info.w, project_info.s, &xmin, &ymin);
- hammer (project_info.e, project_info.n, &xmax, &ymax);
- outside = (PFI) rect_outside;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.check_side = TRUE;
- search = TRUE;
- }
- map_setxy (xmin, xmax, ymin, ymax);
- forward = (PFI)hammer; inverse = (PFI)ihammer;
- gmtdefs.basemap_type = 1;
- return (search);
- }
-
- int vhammer (lon0, scale)
- double lon0, scale; {
- /* Set up Hammer-Aitoff Equal-Area projection */
-
- gmt_check_R_J (&lon0);
- project_info.central_meridian = lon0;
- project_info.w_r = 0.25 * (scale * M_PR_DEG * 360.0); /* = Half the minor axis */
- }
-
- int hammer (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to Hammer-Aitoff Equal-Area x/y */
- double clat, D;
-
- if (fabs(lat) == 90.0) { /* Save time */
- *x = 0.0;
- *y = M_SQRT2 * copysign (EQ_RAD, lat);
- return;
- }
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- lon -= project_info.central_meridian;
- lat *= D2R;
- clat = cos (lat);
- lon *= (0.5 * D2R);
-
- D = EQ_RAD * sqrt (2.0 / (1.0 + clat * cos (lon)));
- *x = 2.0 * D * clat * sin (lon);
- *y = D * sin (lat);
- }
-
- int ihammer (lon, lat, x, y)
- double *lon, *lat, x, y; {
- /* Convert Hammer_Aitoff Equal-Area x/y to lon/lat */
- double rho, c, angle;
-
- x *= 0.5;
- rho = hypot (x, y);
-
- if (rho == 0.0) {
- *lat = 0.0;
- *lon = project_info.central_meridian;
- }
- else {
- c = 2.0 * d_asin (0.5 * rho / EQ_RAD);
- *lat = d_asin (y * sin (c) / rho) * R2D;
- if (fabs (c - M_PI_2) < SMALL)
- angle = (x == 0.0) ? 0.0 : copysign (180.0, x);
- else
- angle = 2.0 * R2D * atan (x * tan (c) / rho);
- *lon = project_info.central_meridian + angle;
- }
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE WINKEL-TRIPEL MODIFIED AZIMUTHAL PROJECTION
- */
-
- int map_init_winkel () {
- int search;
- double xmin, xmax, ymin, ymax, x, y, dummy;
-
- gmt_set_spherical (); /* PW: Force spherical for now */
-
- if (project_info.pars[0] < 0.0) project_info.pars[0] += 360.0;
- gmt_world_map = (fabs (project_info.w - project_info.e) == 360.0);
- if (project_info.gave_map_width) project_info.pars[1] /= (M_PR_DEG * fabs (project_info.w - project_info.e));
- if (project_info.units_pr_degree) project_info.pars[1] /= M_PR_DEG;
- vwinkel (project_info.pars[0], project_info.pars[1]);
- project_info.x_scale = project_info.y_scale = 2.0 * project_info.pars[1] / (1.0 + project_info.r_cosphi1);
-
- if (project_info.region) {
- y = (project_info.s * project_info.n <= 0.0) ? 0.0 : MIN (fabs(project_info.s), fabs(project_info.n));
- x = (fabs (project_info.w - project_info.central_meridian) > fabs (project_info.e - project_info.central_meridian)) ? project_info.w : project_info.e;
- winkel (project_info.w, y, &xmin, &dummy);
- winkel (project_info.e, y, &xmax, &dummy);
- winkel (x, project_info.s, &dummy, &ymin);
- winkel (x, project_info.n, &dummy, &ymax);
- outside = (PFI) wesn_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- left_edge = (PFD) left_winkel;
- right_edge = (PFD) right_winkel;
- frame_info.horizontal = 2;
- search = FALSE;
- }
- else {
- winkel (project_info.w, project_info.s, &xmin, &ymin);
- winkel (project_info.e, project_info.n, &xmax, &ymax);
- outside = (PFI) rect_outside;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.check_side = TRUE;
- search = TRUE;
- }
- map_setxy (xmin, xmax, ymin, ymax);
- forward = (PFI)winkel; inverse = (PFI)iwinkel;
- gmtdefs.basemap_type = 1;
- return (search);
- }
-
- int vwinkel (lon0, scale)
- double lon0, scale; {
- /* Set up Winkel Tripel projection */
-
- gmt_check_R_J (&lon0);
- project_info.r_cosphi1 = cos (50.467 * D2R);
- project_info.central_meridian = lon0;
- project_info.w_r = 0.25 * (scale * M_PR_DEG * 360.0); /* = Half the minor axis */
- }
-
- int winkel (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to Winkel Tripel x/y */
- double C, D, x1, x2, y1, y2;
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- lon -= project_info.central_meridian;
- lat *= D2R;
- lon *= (0.5 * D2R);
-
- /* Fist find Aitoff x/y */
-
- D = d_acos (cos (lat) * cos (lon));
- if (D == 0.0)
- x1 = y1 = 0.0;
- else {
- C = sin (lat) / sin (D);
- D *= EQ_RAD;
- x1 = copysign (2.0 * D * d_sqrt (1.0 - C * C), lon);
- y1 = D * C;
- }
-
- /* Then get equirectangular projection */
-
- x2 = 2.0 * EQ_RAD * lon * project_info.r_cosphi1;
- y2 = EQ_RAD * lat;
-
- /* Winkler is the average value */
-
- *x = 0.5 * (x1 + x2);
- *y = 0.5 * (y1 + y2);
- }
-
- int iwinkel (lon, lat, x, y)
- double *lon, *lat, x, y; {
- /* Convert Winkel Tripel x/y to lon/lat */
- /* Only works if point is on perimeter */
-
- int n_iter = 0;
- double c, phi, phi0, delta;
-
- c = 2.0 * y / EQ_RAD;
- phi = y / EQ_RAD;
- do {
- phi0 = phi;
- phi = phi0 - (phi0 + M_PI_2 * sin (phi0) - c) / (1.0 + M_PI_2 * cos (phi0));
- delta = fabs (phi - phi0);
- n_iter++;
- }
- while (delta > SMALL && n_iter < 100);
- *lat = phi * R2D;
- *lon = project_info.central_meridian + copysign (180.0, x - gmt_half_map_size);
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE ECKERT VI PROJECTION
- */
-
- int map_init_eckert () {
- int search;
- double xmin, xmax, ymin, ymax, y, dummy;
-
- gmt_set_spherical (); /* PW: Force spherical for now */
-
- if (project_info.pars[0] < 0.0) project_info.pars[0] += 360.0;
- gmt_world_map = (fabs (project_info.w - project_info.e) == 360.0);
- if (project_info.gave_map_width) project_info.pars[1] /= (M_PR_DEG * fabs (project_info.w - project_info.e));
- if (project_info.units_pr_degree) project_info.pars[1] /= M_PR_DEG;
- veckert (project_info.pars[0]);
- project_info.x_scale = project_info.y_scale = 0.5 * project_info.pars[1] * sqrt (2.0 + M_PI);
-
- if (project_info.region) {
- y = (project_info.s * project_info.n <= 0.0) ? 0.0 : MIN (fabs(project_info.s), fabs(project_info.n));
- eckert (project_info.w, y, &xmin, &dummy);
- eckert (project_info.e, y, &xmax, &dummy);
- eckert (project_info.central_meridian, project_info.s, &dummy, &ymin);
- eckert (project_info.central_meridian, project_info.n, &dummy, &ymax);
- outside = (PFI) wesn_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- left_edge = (PFD) left_eckert;
- right_edge = (PFD) right_eckert;
- frame_info.horizontal = 2;
- search = FALSE;
- }
- else {
- eckert (project_info.w, project_info.s, &xmin, &ymin);
- eckert (project_info.e, project_info.n, &xmax, &ymax);
- outside = (PFI) rect_outside;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.check_side = TRUE;
- search = TRUE;
- }
- map_setxy (xmin, xmax, ymin, ymax);
- forward = (PFI)eckert; inverse = (PFI)ieckert;
- gmtdefs.basemap_type = 1;
- return (search);
- }
-
- int veckert (lon0)
- double lon0; {
- /* Set up Eckert VI projection */
-
- gmt_check_R_J (&lon0);
- project_info.k_r = EQ_RAD / sqrt (2.0 + M_PI);
- project_info.k_ir = 1.0 / project_info.k_r;
- project_info.central_meridian = lon0;
- }
-
- int eckert (lon, lat, x, y)
- double lon, lat, *x, *y; {
- int n_iter = 0;
- double phi, delta, s_lat;
- /* Convert lon/lat to Eckert VI x/y */
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- lon -= project_info.central_meridian;
- lat *= D2R;
- lon *= D2R;
- phi = lat;
- s_lat = sin (lat);
- do {
- delta = -(phi + sin (phi) - (1.0 + M_PI_2) * s_lat) / (1.0 + cos (phi));
- phi += delta;
- }
- while (delta > SMALL && n_iter < 100);
-
- *x = project_info.k_r * lon * (1.0 + cos (phi));
- *y = 2.0 * project_info.k_r * phi;
- }
-
- int ieckert (lon, lat, x, y)
- double *lon, *lat, x, y; {
- /* Convert Eckert VI x/y to lon/lat */
-
- double phi;
-
- phi = 0.5 * y * project_info.k_ir;
- *lat = d_asin ((phi + sin (phi)) / (1.0 + M_PI_2)) * R2D;
- *lon = project_info.central_meridian + R2D * x * project_info.k_ir / (1.0 + cos (phi));
- }
- /*
- * TRANSFORMATION ROUTINES FOR THE ROBINSON PSEUDOCYLIDRICAL PROJECTION
- */
-
- int map_init_robinson () {
- int search;
- double xmin, xmax, ymin, ymax, y, dummy;
-
- gmt_set_spherical (); /* PW: Force spherical for now */
-
- if (project_info.pars[0] < 0.0) project_info.pars[0] += 360.0;
- gmt_world_map = (fabs (project_info.w - project_info.e) == 360.0);
- if (project_info.gave_map_width) project_info.pars[1] /= (M_PR_DEG * fabs (project_info.w - project_info.e));
- if (project_info.units_pr_degree) project_info.pars[1] /= M_PR_DEG;
- vrobinson (project_info.pars[0]);
- project_info.x_scale = project_info.y_scale = project_info.pars[1] / 0.8487;
-
- if (project_info.region) {
- y = (project_info.s * project_info.n <= 0.0) ? 0.0 : MIN (fabs(project_info.s), fabs(project_info.n));
- robinson (project_info.w, y, &xmin, &dummy);
- robinson (project_info.e, y, &xmax, &dummy);
- robinson (project_info.central_meridian, project_info.s, &dummy, &ymin);
- robinson (project_info.central_meridian, project_info.n, &dummy, &ymax);
- outside = (PFI) wesn_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- left_edge = (PFD) left_robinson;
- right_edge = (PFD) right_robinson;
- frame_info.horizontal = 2;
- search = FALSE;
- }
- else {
- robinson (project_info.w, project_info.s, &xmin, &ymin);
- robinson (project_info.e, project_info.n, &xmax, &ymax);
- outside = (PFI) rect_outside;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.check_side = TRUE;
- search = TRUE;
- }
- map_setxy (xmin, xmax, ymin, ymax);
- forward = (PFI)robinson; inverse = (PFI)irobinson;
- gmtdefs.basemap_type = 1;
- return (search);
- }
-
- int vrobinson (lon0)
- double lon0; {
- /* Set up Robinson projection */
- int i;
- FILE *fp;
- char file[512];
-
- gmt_check_R_J (&lon0);
- project_info.n_cx = 0.8487 * EQ_RAD * D2R;
- project_info.n_cy = 1.3523 * EQ_RAD;
- project_info.central_meridian = lon0;
-
- sprintf (file, "%s/robinson.table\0", LIBDIR);
- if ((fp = fopen (file, "r")) == NULL) {
- fprintf (stderr, "%s: GMT SYNTAX ERROR -Jn -JN OPTION: Could not open interpolation table %s\n", gmt_program, file);
- exit (-1);
- }
- for (i = 0; i < 19; i++) fscanf (fp, "%lf %lf %lf\n", &project_info.n_phi[i], &project_info.n_X[i], &project_info.n_Y[i]);
- fclose (fp);
- }
-
- int robinson (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to Robinson x/y */
- double tmp, X, Y;
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- lon -= project_info.central_meridian;
- tmp = fabs (lat);
-
- intpol (project_info.n_phi, project_info.n_X, 19, 1, &tmp, &X, gmtdefs.interpolant);
- intpol (project_info.n_phi, project_info.n_Y, 19, 1, &tmp, &Y, gmtdefs.interpolant);
- *x = project_info.n_cx * X * lon; /* D2R is in n_cx already */
- *y = project_info.n_cy * copysign (Y, lat);
- }
-
- int irobinson (lon, lat, x, y)
- double *lon, *lat, x, y; {
- /* Convert Robinson x/y to lon/lat */
- double X, Y;
-
- Y = fabs (y / project_info.n_cy);
- intpol (project_info.n_Y, project_info.n_phi, 19, 1, &Y, lat, gmtdefs.interpolant);
- intpol (project_info.n_phi, project_info.n_X, 19, 1, lat, &X, gmtdefs.interpolant);
-
- *lat = copysign (*lat, y);
- *lon = x / (project_info.n_cx * X) + project_info.central_meridian;
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE SINUSOIDAL EQUAL AREA PROJECTION
- */
-
- int map_init_sinusoidal () {
- int search;
- double xmin, xmax, ymin, ymax, dummy, y;
-
- gmt_set_spherical (); /* PW: Force spherical for now */
-
- if (project_info.pars[0] < 0.0) project_info.pars[0] += 360.0;
- gmt_world_map = (fabs (project_info.w - project_info.e) == 360.0);
- if (project_info.gave_map_width) project_info.pars[1] /= (M_PR_DEG * fabs (project_info.w - project_info.e));
- if (project_info.s == -90.0) project_info.edge[0] = FALSE;
- if (project_info.n == 90.0) project_info.edge[2] = FALSE;
- vsinusoidal (project_info.pars[0]);
- if (project_info.units_pr_degree) project_info.pars[1] /= M_PR_DEG;
- project_info.x_scale = project_info.y_scale = project_info.pars[1];
- forward = (PFI)sinusoidal; inverse = (PFI)isinusoidal;
- gmtdefs.basemap_type = 1;
-
- if (project_info.region) {
- y = (project_info.s * project_info.n <= 0.0) ? 0.0 : MIN (fabs(project_info.s), fabs(project_info.n));
- sinusoidal (project_info.central_meridian, project_info.s, &dummy, &ymin);
- sinusoidal (project_info.central_meridian, project_info.n, &dummy, &ymax);
- sinusoidal (project_info.w, y, &xmin, &dummy);
- sinusoidal (project_info.e, y, &xmax, &dummy);
- outside = (PFI) wesn_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- left_edge = (PFD) left_sinusoidal;
- right_edge = (PFD) right_sinusoidal;
- frame_info.horizontal = 2;
- project_info.polar = TRUE;
- search = FALSE;
- }
- else {
- sinusoidal (project_info.w, project_info.s, &xmin, &ymin);
- sinusoidal (project_info.e, project_info.n, &xmax, &ymax);
- outside = (PFI) rect_outside;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.check_side = TRUE;
- search = TRUE;
- }
-
- map_setxy (xmin, xmax, ymin, ymax);
-
- return (search);
- }
-
- int vsinusoidal (lon0)
- double lon0; {
- /* Set up Sinusiodal projection */
-
- gmt_check_R_J (&lon0);
- project_info.central_meridian = lon0;
- }
-
- int sinusoidal (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to Sinusoidal Equal-Area x/y */
-
- lon -= project_info.central_meridian;
- while (lon < -180.0) lon += 360.0;
- while (lon > 180.0) lon -= 360.0;
-
- lon *= D2R;
- lat *= D2R;
-
- *x = EQ_RAD * lon * cos (lat);
- *y = EQ_RAD * lat;
- }
-
- int isinusoidal (lon, lat, x, y)
- double *lon, *lat, x, y; {
- /* Convert Sinusoidal Equal-Area x/y to lon/lat (not implemented) */
- *lat = y / EQ_RAD;
- *lon = ((fabs (fabs (*lat) - M_PI) < SMALL) ? 0.0 : R2D * x / (EQ_RAD * cos (*lat))) + project_info.central_meridian;
- *lat *= R2D;
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE CASSINI PROJECTION
- */
-
- int map_init_cassini () {
- BOOLEAN search;
- double xmin, xmax, ymin, ymax;
-
- vcassini (project_info.pars[0], project_info.pars[1]);
-
- forward = (PFI)cassini; inverse = (PFI)icassini;
- if (project_info.units_pr_degree) project_info.pars[2] /= M_PR_DEG;
- project_info.x_scale = project_info.y_scale = project_info.pars[2];
- gmtdefs.basemap_type = 1;
-
- if (project_info.region) {
- xy_search (&xmin, &xmax, &ymin, &ymax);
- outside = (PFI) wesn_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- left_edge = (PFD) left_albers;
- right_edge = (PFD) right_albers;
- search = FALSE;
- }
- else {
- cassini (project_info.w, project_info.s, &xmin, &ymin);
- cassini (project_info.e, project_info.n, &xmax, &ymax);
- outside = (PFI) rect_outside;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.check_side = TRUE;
- search = TRUE;
- }
-
- frame_info.horizontal = TRUE;
- map_setinfo (xmin, xmax, ymin, ymax, project_info.pars[2]);
-
- return (search);
- }
-
- int vcassini (lon0, lat0)
- double lon0, lat0; {
- /* Set up Cassini projection */
-
- project_info.central_meridian = lon0;
- project_info.pole = lat0;
- lat0 *= D2R;
-
- project_info.c_M0 = EQ_RAD * ((1.0 - 0.25 * ECC2 - 0.046875 * ECC4 - 0.01953125 * ECC6) * lat0
- - (0.375 * ECC2 + 0.09375 * ECC4 + 0.0439453125 * ECC6) * sin (2.0 * lat0)
- + (0.05859375 * ECC4 + 0.0439453125 * ECC6) * sin (4.0 * lat0)
- - 0.0113932292 * ECC6 * sin (6.0 * lat0));
- }
-
- int cassini (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to Cassini x/y */
-
- double tany, cosy, N, T, A, C, M;
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- lon -= project_info.central_meridian;
- lon *= D2R;
- lat *= D2R;
-
- tany = tan (lat);
- cosy = cos (lat);
- N = EQ_RAD / sqrt (1.0 - ECC2 * pow (sin (lat), 2.0));
- T = tany * tany;
- A = lon * cosy;
- C = ECC2 * cosy * cosy / (1.0 - ECC2);
- M = EQ_RAD * ((1.0 - 0.25 * ECC2 - 0.046875 * ECC4 - 0.01953125 * ECC6) * lat
- - (0.375 * ECC2 + 0.09375 * ECC4 + 0.0439453125 * ECC6) * sin (2.0 * lat)
- + (0.05859375 * ECC4 + 0.0439453125 * ECC6) * sin (4.0 * lat)
- - 0.0113932292 * ECC6 * sin (6.0 * lat));
- *x = N * (A - T * pow (A, 3.0) / 6.0 - (8.0 - T + 8 * C) * T * pow (A, 5.0) / 120.0);
- *y = M - project_info.c_M0 + N * tany * (0.5 * A * A + (5.0 - T + 6.0 * C) * pow (A, 4.0) / 24.0);
- }
-
- int icassini (lon, lat, x, y)
- double *lon, *lat, x, y; {
- /* Convert Cassini x/y to lon/lat */
-
- double e1, M1, u1, phi1, tany, siny, T1, N1, R1, D;
-
- e1 = (1.0 - sqrt (1.0 - ECC2)) / (1.0 + sqrt (1.0 - ECC2));
- M1 = project_info.c_M0 + y;
- u1 = M1 / (EQ_RAD * (1.0 - 0.25 * ECC2 - 0.046875 * ECC4 - 0.01953125 * ECC6));
- phi1 = u1 + (1.5 * e1 - 0.84375 * pow (e1, 3.0)) * sin (2.0 * u1)
- + (1.3125 * e1 * e1 - 1.71875 * pow (e1, 4.0)) * sin (4.0 * u1)
- + 1.572916666666 * pow (e1, 3.0) * sin (6.0 * u1)
- + 2.142578125 * pow (e1, 4.0) * sin (8.0 * u1);
- if (fabs (fabs (phi1) - M_PI_2) < SMALL) {
- *lat = copysign (M_PI_2, phi1);
- *lon = project_info.central_meridian;
- }
- else {
- tany = tan (phi1);
- siny = sin (phi1);
- T1 = tany * tany;
- N1 = EQ_RAD / sqrt (1.0 - ECC2 * siny * siny);
- R1 = EQ_RAD * (1.0 - ECC2) / pow (1.0 - ECC2 * siny * siny, 1.5);
- D = x / N1;
-
- *lat = R2D * (phi1 - (N1 * tany / R1) * (0.5 * D * D - (1.0 + 3.0 * T1) * pow (D, 4.0) / 24.0));
- *lon = project_info.central_meridian + R2D * (D - T1 * pow (D, 3.0) / 3.0 + (1.0 + 3.0 * T1) * T1 * pow (D, 5.0)/ 15.0) / cos (phi1);
- }
- }
-
- /*
- * TRANSFORMATION ROUTINES FOR THE ALBERS PROJECTION
- */
-
- int map_init_albers () {
- BOOLEAN search;
- double xmin, xmax, ymin, ymax;
-
- valbers (project_info.pars[0], project_info.pars[1], project_info.pars[2], project_info.pars[3]);
- if (project_info.units_pr_degree) project_info.pars[4] /= M_PR_DEG;
- forward = (PFI)albers; inverse = (PFI)ialbers;
- project_info.x_scale = project_info.y_scale = project_info.pars[4];
-
- if (project_info.region) {
- xy_search (&xmin, &xmax, &ymin, &ymax);
- outside = (PFI) wesn_outside;
- crossing = (PFI) wesn_crossing;
- overlap = (PFI) wesn_overlap;
- map_clip = (PFI) wesn_clip;
- left_edge = (PFD) left_albers;
- right_edge = (PFD) right_albers;
- search = FALSE;
- }
- else {
- albers (project_info.w, project_info.s, &xmin, &ymin);
- albers (project_info.e, project_info.n, &xmax, &ymax);
- outside = (PFI) rect_outside;
- crossing = (PFI) rect_crossing;
- overlap = (PFI) rect_overlap;
- map_clip = (PFI) rect_clip;
- left_edge = (PFD) left_rect;
- right_edge = (PFD) right_rect;
- frame_info.check_side = TRUE;
- search = TRUE;
- }
- frame_info.horizontal = TRUE;
- gmtdefs.n_lat_nodes = 2;
- map_setinfo (xmin, xmax, ymin, ymax, project_info.pars[4]);
-
- return (search);
- }
-
-
- int valbers (lon0, lat0, ph1, ph2)
- double lon0, lat0, ph1, ph2; {
- /* Set up Albers projection */
-
- double s0, s1, s2, q0, q1, q2, m1, m2;
-
- project_info.central_meridian = lon0;
- project_info.north_pole = (lat0 > 0.0);
- project_info.pole = (project_info.north_pole) ? 90.0 : -90.0;
- lat0 *= D2R;
- ph1 *= D2R; ph2 *= D2R;
-
- s0 = sin (lat0); s1 = sin (ph1); s2 = sin (ph2);
- m1 = cos (ph1) * cos (ph1) / (1.0 - ECC2 * s1 * s1); /* Actually m1 and m2 squared */
- m2 = cos (ph2) * cos (ph2) / (1.0 - ECC2 * s2 * s2);
- q0 = (1.0 - ECC2) * (s0 / (1.0 - ECC2 * s0 * s0) - (0.5 / ECC) * log ((1.0 - ECC * s0) / (1.0 + ECC * s0)));
- q1 = (1.0 - ECC2) * (s1 / (1.0 - ECC2 * s1 * s1) - (0.5 / ECC) * log ((1.0 - ECC * s1) / (1.0 + ECC * s1)));
- q2 = (1.0 - ECC2) * (s2 / (1.0 - ECC2 * s2 * s2) - (0.5 / ECC) * log ((1.0 - ECC * s2) / (1.0 + ECC * s2)));
-
- project_info.a_n = (m1 - m2) / (q2 - q1);
- project_info.a_C = m1 + project_info.a_n * q1;
- project_info.a_rho0 = EQ_RAD * sqrt (project_info.a_C - project_info.a_n * q0) / project_info.a_n;
- }
-
- int albers (lon, lat, x, y)
- double lon, lat, *x, *y; {
- /* Convert lon/lat to Albers x/y */
-
- double s, q, theta, rho;
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- lon -= project_info.central_meridian;
- lon *= D2R;
- lat *= D2R;
-
- s = sin (lat);
- q = (1.0 - ECC2) * (s / (1.0 - ECC2 * s * s) - (0.5 / ECC) * log ((1.0 - ECC * s) / (1.0 + ECC * s)));
- theta = project_info.a_n * lon;
- rho = EQ_RAD * sqrt (project_info.a_C - project_info.a_n * q) / project_info.a_n;
-
- *x = rho * sin (theta);
- *y = project_info.a_rho0 - rho * cos (theta);
- }
-
- int ialbers (lon, lat, x, y)
- double *lon, *lat, x, y; {
- /* Convert ALbers x/y to lon/lat */
-
- int n_iter = 0;
- double theta, rho, q, phi, phi0, s, delta, test;
-
- theta = (project_info.a_n < 0.0) ? d_atan2 (-x, y - project_info.a_rho0) : d_atan2 (x, project_info.a_rho0 - y);
- rho = hypot (x, project_info.a_rho0 - y);
- q = (project_info.a_C - rho * rho * project_info.a_n * project_info.a_n / (EQ_RAD * EQ_RAD)) / project_info.a_n;
-
- test = 1.0 - (0.5 * (1.0 - ECC2) / ECC) * log ((1.0 - ECC) / (1.0 + ECC));
- if (fabs (fabs (q) - test) < SMALL)
- *lat = copysign (90.0, q);
- else {
- phi = d_asin (0.5 * q);
- do {
- phi0 = phi;
- s = sin (phi0);
- phi = phi0 + 0.5 * (1.0 - ECC2 * s * s) * ((q / (1.0 - ECC2)) - s / (1.0 - ECC2 * s * s)
- + 0.5 * log ((1 - ECC * s) / (1.0 + ECC * s)) / ECC) / cos (phi);
- delta = fabs (phi - phi0);
- n_iter++;
- }
- while (delta > SMALL && n_iter < 100);
- *lat = R2D * phi;
- }
- *lon = project_info.central_meridian + R2D * theta / project_info.a_n;
- }
-
- /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- *
- * S E C T I O N 1.1 : S U P P O R T I N G R O U T I N E S
- *
- * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
-
- int wesn_search () {
- double dx, dy, w, e, s, n, x, y, lon, lat;
- int i, j;
-
- /* Search for extreme original coordinates lon/lat */
-
- dx = project_info.xmax / gmtdefs.n_lon_nodes;
- dy = project_info.ymax / gmtdefs.n_lat_nodes;
- w = s = MAXDOUBLE; e = n = -MAXDOUBLE;
- for (i = 0; i <= gmtdefs.n_lon_nodes; i++) {
- x = i * dx;
- xy_to_geo (&lon, &lat, x, 0.0);
- if (lon < w) w = lon;
- if (lon > e) e = lon;
- if (lat < s) s = lat;
- if (lat > n) n = lat;
- xy_to_geo (&lon, &lat, x, project_info.ymax);
- if (lon < w) w = lon;
- if (lon > e) e = lon;
- if (lat < s) s = lat;
- if (lat > n) n = lat;
- }
- for (j = 0; j <= gmtdefs.n_lat_nodes; j++) {
- y = j * dy;
- xy_to_geo (&lon, &lat, 0.0, y);
- if (lon < w) w = lon;
- if (lon > e) e = lon;
- if (lat < s) s = lat;
- if (lat > n) n = lat;
- xy_to_geo (&lon, &lat, project_info.xmax, y);
- if (lon < w) w = lon;
- if (lon > e) e = lon;
- if (lat < s) s = lat;
- if (lat > n) n = lat;
- }
-
- /* Then check if one or both poles are inside map; then the above wont be correct */
-
- if (!map_outside (project_info.central_meridian, 90.0)) {
- n = 90.0;
- w = 0.0;
- e = 360.0;
- }
- /* if (project_info.projection != AZ_EQDIST && !map_outside (project_info.central_meridian, -90.0)) { why was it like this? */
- if (!map_outside (project_info.central_meridian, -90.0)) {
- s = -90.0;
- w = 0.0;
- e = 360.0;
- }
-
- s -= 0.1; if (s < -90.0) s = -90.0; /* Make sure point is not inside area, 0.1 is just a small number */
- n += 0.1; if (n > 90.0) n = 90.0;
- w -= 0.1; e += 0.1;
-
- if (fabs (w - e) > 360.0) {
- w = 0.0;
- e = 360.0;
- }
- project_info.w = w;
- project_info.e = e;
- project_info.s = s;
- project_info.n = n;
- }
-
- int xy_search (x0, x1, y0, y1)
- double *x0, *x1, *y0, *y1; {
- int i, j;
- double xmin, xmax, ymin, ymax, w, s, x, y;
-
- /* Find min/max forward values */
-
- xmax = ymax = -MAXDOUBLE;
- xmin = ymin = MAXDOUBLE;
- gmtdefs.dlon = fabs (project_info.e - project_info.w) / 500;
- gmtdefs.dlat = fabs (project_info.n - project_info.s) / 500;
-
- for (i = 0; i <= 500; i++) {
- w = project_info.w + i * gmtdefs.dlon;
- (*forward) (w, project_info.s, &x, &y);
- if (x < xmin) xmin = x;
- if (y < ymin) ymin = y;
- if (x > xmax) xmax = x;
- if (y > ymax) ymax = y;
- (*forward) (w, project_info.n, &x, &y);
- if (x < xmin) xmin = x;
- if (y < ymin) ymin = y;
- if (x > xmax) xmax = x;
- if (y > ymax) ymax = y;
- }
- for (j = 0; j <= 500; j++) {
- s = project_info.s + j * gmtdefs.dlat;
- (*forward) (project_info.w, s, &x, &y);
- if (x < xmin) xmin = x;
- if (y < ymin) ymin = y;
- if (x > xmax) xmax = x;
- if (y > ymax) ymax = y;
- (*forward) (project_info.e, s, &x, &y);
- if (x < xmin) xmin = x;
- if (y < ymin) ymin = y;
- if (x > xmax) xmax = x;
- if (y > ymax) ymax = y;
- }
-
- *x0 = xmin; *x1 = xmax; *y0 = ymin; *y1 = ymax;
- }
-
- int map_crossing (lon1, lat1, lon2, lat2, xlon, xlat, xx, yy, sides)
- double lon1, lat1, lon2, lat2, *xlon, *xlat, *xx, *yy;
- int *sides; {
- int nx;
- gmt_corner = -1;
- nx = (*crossing) (lon1, lat1, lon2, lat2, xlon, xlat, xx, yy, sides);
-
- /* nx may be -2, in which case we dont want to check the order */
- if (nx == 2) { /* Must see if crossings are in correct order */
- double da, db;
-
- da = great_circle_dist (lon1, lat1, xlon[0], xlat[0]);
- db = great_circle_dist (lon1, lat1, xlon[1], xlat[1]);
- if (da > db) { /* Must swap */
- d_swap (xlon[0], xlon[1]);
- d_swap (xlat[0], xlat[1]);
- d_swap (xx[0], xx[1]);
- d_swap (yy[0], yy[1]);
- i_swap (sides[0], sides[1]);
- }
- }
- return (abs(nx));
- }
-
- int map_outside (lon, lat)
- double lon, lat; {
- gmt_x_status_old = gmt_x_status_new;
- gmt_y_status_old = gmt_y_status_new;
- return ((*outside) (lon, lat));
- }
-
- int break_through (x0, y0, x1, y1)
- double x0, y0, x1, y1; {
-
- if (gmt_x_status_old == gmt_x_status_new && gmt_y_status_old == gmt_y_status_new) return (FALSE);
- if ((gmt_x_status_old == 0 && gmt_y_status_old == 0) || (gmt_x_status_new == 0 && gmt_y_status_new == 0)) return (TRUE);
-
- /* Less clearcut case, check for overlap */
-
- return ( (*overlap) (x0, y0, x1, y1));
- }
-
- int rect_overlap (lon0, lat0, lon1, lat1)
- double lon0, lat0, lon1, lat1; {
- double x0, y0, x1, y1;
-
- geo_to_xy (lon0, lat0, &x0, &y0);
- geo_to_xy (lon1, lat1, &x1, &y1);
-
- if (x0 > x1) d_swap (x0, x1);
- if (y0 > y1) d_swap (y0, y1);
- if (x1 < project_info.xmin || x0 > project_info.xmax) return (FALSE);
- if (y1 < project_info.ymin || y0 > project_info.ymax) return (FALSE);
- return (TRUE);
- }
-
- int wesn_overlap (lon0, lat0, lon1, lat1)
- double lon0, lat0, lon1, lat1; {
- if (lon0 > lon1) d_swap (lon0, lon1);
- if (lat0 > lat1) d_swap (lat0, lat1);
- if (lon1 < project_info.w) {
- lon0 += 360.0;
- lon1 += 360.0;
- }
- else if (lon0 > project_info.e) {
- lon0 -= 360.0;
- lon1 -= 360.0;
- }
- if (lon1 < project_info.w || lon0 > project_info.e) return (FALSE);
- if (lat1 < project_info.s || lat0 > project_info.n) return (FALSE);
- return (TRUE);
- }
-
- int radial_overlap (lon0, lat0, lon1, lat1)
- double lon0, lat0, lon1, lat1; { /* Dummy routine */
- return (TRUE);
- }
-
- int wrap_around_check (angle, last_x, last_y, this_x, this_y, xx, yy, sides, nx)
- double *angle, last_x, last_y, this_x, this_y, *xx, *yy;
- int *sides, *nx; {
- int i, wrap = FALSE, skip;
- double dx, width, jump, gmt_half_map_width();
-
- jump = this_x - last_x;
- width = gmt_half_map_width (this_y);
-
- skip = (fabs (jump) < width || (fabs(jump) <= SMALL || fabs(width) <= SMALL));
-
- for (i = 0; i < (*nx); i++) { /* Must check if the crossover found should wrap around */
- if (skip) continue;
-
- if (jump < (-width)) { /* Crossed right boundary */
- xx[0] = right_boundary (this_y); xx[1] = left_boundary (this_y);
- dx = this_x + gmt_map_width - last_x;
- yy[0] = yy[1] = last_y + (gmt_map_width - last_x) * (this_y - last_y) / dx;
- sides[0] = 1; sides[1] = 3;
- angle[0] = R2D * d_atan2 (this_y - last_y, dx);
- }
- else { /* Crossed left boundary */
- xx[0] = left_boundary (this_y); xx[1] = right_boundary (this_y);
- dx = last_x + gmt_map_width - this_x;
- yy[0] = yy[1] = last_y + last_x * (this_y - last_y) / dx;
- sides[0] = 3; sides[1] = 1;
- angle[0] = R2D * d_atan2 (this_y - last_y, -dx);
- }
- angle[1] = angle[0] + 180.0;
- if (yy[0] >= 0.0 && yy[0] <= project_info.ymax) wrap = TRUE;
- }
-
- if (*nx == 0 && !skip) { /* Must wrap around */
- if (jump < (-width)) { /* Crossed right boundary */
- xx[0] = right_boundary (this_y); xx[1] = left_boundary (this_y);
- dx = this_x + gmt_map_width - last_x;
- yy[0] = yy[1] = last_y + (gmt_map_width - last_x) * (this_y - last_y) / dx;
- sides[0] = 1; sides[1] = 3;
- angle[0] = R2D * d_atan2 (this_y - last_y, dx);
- }
- else { /* Crossed left boundary */
- xx[0] = left_boundary (this_y); xx[1] = right_boundary (this_y);
- dx = last_x + gmt_map_width - this_x;
- yy[0] = yy[1] = last_y + last_x * (this_y - last_y) / dx;
- sides[0] = 3; sides[1] = 1;
- angle[0] = R2D * d_atan2 (this_y - last_y, -dx);
- }
- if (yy[0] >= 0.0 && yy[0] <= project_info.ymax) wrap = TRUE;
- angle[1] = angle[0] + 180.0;
- }
-
- if (wrap) *nx = 2;
- return (wrap);
- }
-
- double left_ellipse (y)
- double y; {
- y -= project_info.w_r;
- return (gmt_half_map_size - 2.0 * d_sqrt (project_info.w_r * project_info.w_r - y * y));
- }
-
- double left_winkel (y)
- double y; {
- int n_iter = 0;
- double c, phi, phi0, delta, x, y0;
-
- y0 = 0.5 * project_info.ymax;
- y -= y0;
- y /= project_info.y_scale;
- c = 2.0 * y / EQ_RAD;
- phi = y / EQ_RAD;
- do {
- phi0 = phi;
- phi = phi0 - (phi0 + M_PI_2 * sin (phi0) - c) / (1.0 + M_PI_2 * cos (phi0));
- delta = fabs (phi - phi0);
- n_iter++;
- }
- while (delta > SMALL && n_iter < 100);
- geo_to_xy (project_info.central_meridian-180.0, phi * R2D, &x, &c);
- return (x);
- }
-
- double left_eckert (y)
- double y; {
- double x, phi;
-
- y -= project_info.y0;
- y /= project_info.y_scale;
- phi = 0.5 * y * project_info.k_ir;
- x = project_info.k_r * D2R * (project_info.w - project_info.central_meridian) * (1.0 + cos (phi));
- return (x * project_info.x_scale + project_info.x0);
- }
-
- double left_robinson (y)
- double y; {
- double x, X, Y;
-
- y -= project_info.y0;
- y /= project_info.y_scale;
- Y = fabs (y / project_info.n_cy);
- intpol (project_info.n_Y, project_info.n_X, 19, 1, &Y, &X, gmtdefs.interpolant);
-
- x = project_info.n_cx * X * (project_info.w - project_info.central_meridian);
- return (x * project_info.x_scale + project_info.x0);
- }
-
- double left_sinusoidal (y)
- double y; {
- double x, lat;
- y -= project_info.y0;
- y /= project_info.y_scale;
- lat = y / EQ_RAD;
- x = EQ_RAD * (project_info.w - project_info.central_meridian) * D2R * cos (lat);
- return (x * project_info.x_scale + project_info.x0);
- }
-
- double left_circle (y)
- double y; {
- y -= project_info.r;
- return (gmt_half_map_size - d_sqrt (project_info.r * project_info.r - y * y));
- }
-
- double left_lambert (y)
- double y; {
- fprintf (stderr, "GMT Warning: left_lambert () called which should not occur. Report to gmt@soest\n");
- return (0.0); /* For completeness, but should never be called */
- }
-
- double left_albers (y)
- double y; {
- fprintf (stderr, "GMT Warning: left_albers () called which should not occur. Report to gmt@soest\n");
- return (0.0); /* For completeness, but should never be called */
- }
-
- double left_rect (y)
- double y; {
- return (0.0);
- }
-
- double left_boundary (y)
- double y; {
- return ((*left_edge) (y));
- }
-
- double right_ellipse (y)
- double y; {
- y -= project_info.w_r;
- return (gmt_half_map_size + 2.0 * d_sqrt (project_info.w_r * project_info.w_r - y * y));
- }
-
- double right_winkel (y)
- double y; {
- int n_iter = 0;
- double c, phi, phi0, delta, x, y0;
-
- y0 = 0.5 * project_info.ymax;
- y -= y0;
- y /= project_info.y_scale;
- c = 2.0 * y / EQ_RAD;
- phi = y / EQ_RAD;
- do {
- phi0 = phi;
- phi = phi0 - (phi0 + M_PI_2 * sin (phi0) - c) / (1.0 + M_PI_2 * cos (phi0));
- delta = fabs (phi - phi0);
- n_iter++;
- }
- while (delta > SMALL && n_iter < 100);
- geo_to_xy (project_info.central_meridian+180.0, phi * R2D, &x, &c);
- return (x);
- }
-
- double right_eckert (y)
- double y; {
- double x, phi;
-
- y -= project_info.y0;
- y /= project_info.y_scale;
- phi = 0.5 * y * project_info.k_ir;
- x = project_info.k_r * D2R * (project_info.e - project_info.central_meridian) * (1.0 + cos (phi));
- return (x * project_info.x_scale + project_info.x0);
- }
-
- double right_robinson (y)
- double y; {
- double x, X, Y;
-
- y -= project_info.y0;
- y /= project_info.y_scale;
- Y = fabs (y / project_info.n_cy);
- intpol (project_info.n_Y, project_info.n_X, 19, 1, &Y, &X, gmtdefs.interpolant);
-
- x = project_info.n_cx * X * (project_info.e - project_info.central_meridian);
- return (x * project_info.x_scale + project_info.x0);
- }
-
- double right_sinusoidal (y)
- double y; {
- double x, lat;
- y -= project_info.y0;
- y /= project_info.y_scale;
- lat = y / EQ_RAD;
- x = EQ_RAD * (project_info.e - project_info.central_meridian) * D2R * cos (lat);
- return (x * project_info.x_scale + project_info.x0);
- }
-
- double right_circle (y)
- double y; {
- y -= project_info.r;
- return (gmt_half_map_size + d_sqrt (project_info.r * project_info.r - y * y));
- }
-
- double right_lambert (y)
- double y; {
- fprintf (stderr, "GMT Warning: right_lambert () called which should not occur. Report to gmt@soest\n");
- return (gmt_map_width); /* For completeness, but should never be called */
- }
-
- double right_albers (y)
- double y; {
- fprintf (stderr, "GMT Warning: right_albers () called which should not occur. Report to gmt@soest\n");
- return (gmt_map_width); /* For completeness, but should never be called */
- }
-
- double right_rect (y)
- double y; {
- return (gmt_map_width);
- }
-
- double right_boundary (y)
- double y; {
- return ((*right_edge) (y));
- }
-
- int map_jump (x0, y0, x1, y1)
- double x0, y0, x1, y1; {
- /* TRUE if x-distance between points exceeds 1/2 map width at this y value */
- double dx, map_half_size;
-
- if (!MAPPING || fabs (project_info.w - project_info.e) < 90.0) return (FALSE);
-
- map_half_size = gmt_half_map_width (y0);
- if (fabs (map_half_size) < SMALL) return (0);
-
- dx = x1 - x0;
- if (dx > map_half_size) /* Cross left/west boundary */
- return (-1);
- else if (dx < (-map_half_size)) /* Cross right/east boundary */
- return (1);
- else
- return (0);
- }
-
- int get_crossings (xc, yc, x0, y0, x1, y1)
- double xc[], yc[];
- double x0, y0, x1, y1; { /* Finds crossings for wrap-arounds */
- double xa, xb, ya, yb, dxa, dxb, dyb, c;
-
- xa = x0; xb = x1;
- ya = y0; yb = y1;
- if (xa > xb) {
- d_swap (xa, xb);
- d_swap (ya, yb);
- }
-
- xb -= 2.0 * gmt_half_map_width (yb);
-
- dxa = xa - left_boundary (ya);
- dxb = left_boundary (yb) - xb;
- c = (dxb == 0.0) ? 0.0 : 1.0 + dxa/dxb;
- dyb = (c == 0.0) ? 0.0 : fabs (yb - ya) / c;
- yc[0] = yc[1] = (ya > yb) ? yb + dyb : yb - dyb;
- xc[0] = left_boundary (yc[0]);
- xc[1] = right_boundary (yc[0]);
- }
-
- double gmt_half_map_width (y)
- double y; { /* Returns 1/2-width of map in inches given y value */
- double half_width;
-
- switch (project_info.projection) {
-
- case STEREO: /* Must compute width of circular map based on y value */
- case LAMB_AZ_EQ:
- case ORTHO:
- case AZ_EQDIST:
- if (project_info.region && gmt_world_map) {
- y -= project_info.r;
- half_width = d_sqrt (project_info.r * project_info.r - y * y);
- }
- else
- half_width = gmt_half_map_size;
- break;
-
- case MOLLWEIDE: /* Must compute width of Mollweide map based on y value */
- case HAMMER:
- if (project_info.region && gmt_world_map) {
- y -= project_info.w_r;
- half_width = 2.0 * d_sqrt (project_info.w_r * project_info.w_r - y * y);
- }
- else
- half_width = gmt_half_map_size;
- break;
- case WINKEL:
- case SINUSOIDAL:
- case ROBINSON:
- case ECKERT:
- if (project_info.region && gmt_world_map)
- half_width = right_boundary (y) - gmt_half_map_size;
- else
- half_width = gmt_half_map_size;
- break;
-
- default: /* Rectangular maps are easy */
- half_width = gmt_half_map_size;
- break;
- }
- return (half_width);
- }
-
- int will_it_wrap (x, y, n, start)
- double x[], y[];
- int n, *start; { /* Determines if a polygon will wrap at edges */
- int i;
- BOOLEAN wrap;
- double dx, w;
-
- if (!gmt_world_map) return (FALSE);
-
- for (i = 1, wrap = FALSE; !wrap && i < n; i++) {
- dx = fabs (x[i] - x[i-1]);
- if (dx > (w = gmt_half_map_width (y[i])) && w > 0.0)
- wrap = TRUE;
- }
- *start = i - 1;
- return (wrap);
- }
-
- int truncate_left (x, y, n, start)
- double x[], y[];
- int n, start; { /* Truncates a wrapping polygon agains left edge */
- int i, i1, j, k;
- BOOLEAN on;
- double dx, xc[2], yc[2];
-
- /* Find first point that is left of map center */
-
- i = (x[start] < gmt_half_map_size) ? start : start - 1;
-
- if (!gmt_n_alloc) get_plot_array ();
-
- gmt_x_plot[0] = x[i]; gmt_y_plot[0] = y[i];
- for (k = j = 1, on = FALSE; k < n; k++, j++) {
- i1 = i;
- i = (i + 1)%n; /* Next point */
- dx = fabs (x[i] - x[i1]);
- if (dx > gmt_half_map_width (y[i])) {
- on = !on;
- get_crossings (xc, yc, x[i1], y[i1], x[i], y[i]);
- gmt_x_plot[j] = xc[0];
- gmt_y_plot[j] = yc[0];
- j++;
- if (j >= gmt_n_alloc) get_plot_array ();
- }
- gmt_x_plot[j] = (on) ? left_boundary (y[i]) : x[i];
- gmt_y_plot[j] = y[i];
- if (j >= gmt_n_alloc) get_plot_array ();
- }
- return (j);
- }
-
- int truncate_right (x, y, n, start)
- double x[], y[];
- int n, start; { /* Truncates a wrapping polygon agains right edge */
- int i, i1, j, k;
- BOOLEAN on;
- double dx, xc[2], yc[2];
-
- /* Find first point that is right of map center */
-
- i = (x[start] > gmt_half_map_size) ? start : start - 1;
-
- if (!gmt_n_alloc) get_plot_array ();
- gmt_x_plot[0] = x[i]; gmt_y_plot[0] = y[i];
- for (k = j = 1, on = FALSE; k < n; k++, j++) {
- i1 = i;
- i = (i + 1)%n; /* Next point */
- dx = fabs (x[i] - x[i1]);
- if (dx > gmt_half_map_width (y[i])) {
- on = !on;
- get_crossings (xc, yc, x[i1], y[i1], x[i], y[i]);
- gmt_x_plot[j] = xc[1];
- gmt_y_plot[j] = yc[1];
- j++;
- if (j >= gmt_n_alloc) get_plot_array ();
- }
- gmt_x_plot[j] = (on) ? right_boundary (y[i]) : x[i];
- gmt_y_plot[j] = y[i];
- if (j >= gmt_n_alloc) get_plot_array ();
- }
- return (j);
- }
-
- double great_circle_dist( lon1, lat1, lon2, lat2 )
- double lon1, lat1, lon2, lat2;
- {
- /* great circle distance on a sphere in degrees */
- double C, a, b, c;
- double cosC, cosa, cosb, cosc;
- double sina, sinb;
-
- if( (lat1==lat2) && (lon1==lon2) ) {
- return( (double)0.0 );
- }
-
- a=D2R*(90.0-lat2);
- b=D2R*(90.0-lat1);
-
- C = D2R*( lon2 - lon1 );
-
- cosa = cos(a);
- cosb = cos(b);
- sina = sin(a);
- sinb = sin(b);
- cosC = cos(C);
-
- cosc = cosa*cosb + sina*sinb*cosC;
- if (cosc<-1.0) c=M_PI; else if (cosc>1) c=0.0; else c=d_acos(cosc);
-
- return( c * R2D);
- }
-
- /* The *_outside rutines returns the status of the current point. Status is
- * the sum of x_status and y_status. x_status may be
- * 0 w < lon < e
- * -1 lon == w
- * 1 lon == e
- * -2 lon < w
- * 2 lon > e
- * y_status may be
- * 0 s < lat < n
- * -1 lat == s
- * 1 lat == n
- * -2 lat < s
- * 2 lat > n
- */
-
- int wesn_outside (lon, lat)
- double lon, lat; {
-
- if (gmt_world_map) {
- while (lon < project_info.w) lon += 360.0;
- while (lon > project_info.e) lon -= 360.0;
- }
-
- if (on_border_is_outside && fabs (lon - project_info.w) < SMALL )
- gmt_x_status_new = -1;
- else if (on_border_is_outside && fabs (lon - project_info.e) < SMALL)
- gmt_x_status_new = 1;
- else if (lon < project_info.w)
- gmt_x_status_new = -2;
- else if (lon > project_info.e)
- gmt_x_status_new = 2;
- else
- gmt_x_status_new = 0;
-
- if (on_border_is_outside && fabs (lat - project_info.s) < SMALL )
- gmt_y_status_new = -1;
- else if (on_border_is_outside && fabs (lat - project_info.n) < SMALL)
- gmt_y_status_new = 1;
- else if (lat < project_info.s )
- gmt_y_status_new = -2;
- else if (lat > project_info.n)
- gmt_y_status_new = 2;
- else
- gmt_y_status_new = 0;
-
- return ( !(gmt_x_status_new == 0 && gmt_y_status_new == 0));
-
- }
-
- int polar_outside (lon, lat)
- double lon, lat; {
-
- if (on_border_is_outside && fabs (lon - project_info.w) < SMALL )
- gmt_x_status_new = -1;
- else if (on_border_is_outside && fabs (lon - project_info.e) < SMALL)
- gmt_x_status_new = 1;
- else if (lon < project_info.w)
- gmt_x_status_new = -2;
- else if (lon > project_info.e)
- gmt_x_status_new = 2;
- else
- gmt_x_status_new = 0;
- if (!project_info.edge[1]) gmt_x_status_new = 0; /* 360 degrees, no edge */
-
- if (on_border_is_outside && fabs (lat - project_info.s) < SMALL )
- gmt_y_status_new = -1;
- else if (on_border_is_outside && fabs (lat - project_info.n) < SMALL)
- gmt_y_status_new = 1;
- else if (lat < project_info.s )
- gmt_y_status_new = -2;
- else if (lat > project_info.n)
- gmt_y_status_new = 2;
- else
- gmt_y_status_new = 0;
- if (gmt_y_status_new < 0 && !project_info.edge[0]) gmt_y_status_new = 0; /* South pole enclosed */
- if (gmt_y_status_new > 0 && !project_info.edge[2]) gmt_y_status_new = 0; /* North pole enclosed */
-
- return ( !(gmt_x_status_new == 0 && gmt_y_status_new == 0));
-
- }
-
- int eqdist_outside (lon, lat)
- double lon, lat; {
- double dlon, cc;
-
- while ((lon - project_info.central_meridian) < -180.0) lon += 360.0;
- while ((lon - project_info.central_meridian) > 180.0) lon -= 360.0;
- dlon = (lon - project_info.central_meridian);
- cc = project_info.sinp * sind (lat) + project_info.cosp * cosd (lat) * cosd (dlon);
- /* if (cc <= -1.0) { */
- if (cc < -1.0) {
- gmt_y_status_new = -1;
- gmt_x_status_new = 0;
- }
- else
- gmt_x_status_new = gmt_y_status_new = 0;
- return ( !(gmt_y_status_new == 0));
- }
-
- int radial_outside (lon, lat)
- double lon, lat; {
- double dist;
-
- /* Test if point is more than 90 spherical degrees from origin. For global maps, let all borders by "south" */
-
- gmt_x_status_new = 0;
- dist = great_circle_dist (lon, lat, project_info.central_meridian, project_info.pole);
- if (on_border_is_outside && fabs (dist - 90.0) < SMALL)
- gmt_y_status_new = -1;
- else if (dist > 90.0)
- gmt_y_status_new = -2;
- else
- gmt_y_status_new = 0;
- return ( !(gmt_y_status_new == 0));
- }
-
- int rect_outside (lon, lat)
- double lon, lat; {
- double x, y;
-
- geo_to_xy (lon, lat, &x, &y);
-
- if (on_border_is_outside && fabs (x - project_info.xmin) < SMALL )
- gmt_x_status_new = -1;
- else if (on_border_is_outside && fabs (x - project_info.xmax) < SMALL)
- gmt_x_status_new = 1;
- else if (x < project_info.xmin)
- gmt_x_status_new = -2;
- else if (x > project_info.xmax)
- gmt_x_status_new = 2;
- else
- gmt_x_status_new = 0;
-
- if (on_border_is_outside && fabs (y -project_info.ymin) < SMALL )
- gmt_y_status_new = -1;
- else if (on_border_is_outside && fabs (y - project_info.ymax) < SMALL)
- gmt_y_status_new = 1;
- else if (y < project_info.ymin)
- gmt_y_status_new = -2;
- else if (y > project_info.ymax)
- gmt_y_status_new = 2;
- else
- gmt_y_status_new = 0;
-
- return ( !(gmt_x_status_new == 0 && gmt_y_status_new == 0));
-
- }
-
- int rect_outside2 (lon, lat)
- double lon, lat; { /* For Azimuthal proj with rect borders since rect_outside may fail for antipodal points */
- if (radial_outside (lon, lat)) return (TRUE); /* Point > 90 degrees away */
- return (rect_outside (lon, lat)); /* Must check if inside box */
- }
-
- int pen_status () {
- int pen = 3;
-
- if (gmt_x_status_old == 0 && gmt_y_status_old == 0)
- pen = 2;
- else if (gmt_x_status_new == 0 && gmt_y_status_new == 0)
- pen = 3;
- return (pen);
- }
-
- int wesn_crossing (lon0, lat0, lon1, lat1, clon, clat, xx, yy, sides)
- double lon0, lat0, lon1, lat1, *clon, *clat, *xx, *yy;
- int *sides; {
- /* Compute the crossover point(s) on the map boundary for rectangular projections */
- int n = 0, i;
- double dlat, dlon;
-
- /* Since it may not be obvious which side the line may cross, and since in some cases the two points may be
- * entirely outside the region but still cut through it, we first find all possible candidates and then decide
- * which ones are valid crossings. We may find 0, 1, or 2 intersections */
-
- /* First align the longitudes to the region */
-
- if (gmt_world_map) {
- while (lon0 < project_info.w) lon0 += 360.0;
- while (lon0 > project_info.e) lon0 -= 360.0;
- while (lon1 < project_info.w) lon1 += 360.0;
- while (lon1 > project_info.e) lon1 -= 360.0;
- }
-
- /* Then set 'almost'-corners to corners */
-
- x_wesn_corner (&lon0);
- x_wesn_corner (&lon1);
- y_wesn_corner (&lat0);
- y_wesn_corner (&lat1);
-
-
- /* if ((lat0 > project_info.s && lat1 <= project_info.s) || (lat1 > project_info.s && lat0 <= project_info.s)) { */
- if ((lat0 >= project_info.s && lat1 <= project_info.s) || (lat1 >= project_info.s && lat0 <= project_info.s)) {
- sides[n] = 0;
- clat[n] = project_info.s;
- dlat = lat0 - lat1;
- clon[n] = (dlat == 0.0) ? lon1 : lon1 + (lon0 - lon1) * (clat[n] - lat1) / dlat;
- x_wesn_corner (&clon[n]);
- if (fabs(dlat) > 0.0 && lon_inside (clon[n], project_info.w, project_info.e)) n++;
- }
- /* if ((lon0 >= project_info.e && lon1 < project_info.e) || (lon1 >= project_info.e && lon0 < project_info.e)) { */
- if ((lon0 >= project_info.e && lon1 <= project_info.e) || (lon1 >= project_info.e && lon0 <= project_info.e)) {
- sides[n] = 1;
- clon[n] = project_info.e;
- dlon = lon0 - lon1;
- clat[n] = (dlon == 0.0) ? lat1 : lat1 + (lat0 - lat1) * (clon[n] - lon1) / dlon;
- y_wesn_corner (&clat[n]);
- if (fabs(dlon) > 0.0 && clat[n] >= project_info.s && clat[n] <= project_info.n) n++;
- }
- /* if ((lat0 >= project_info.n && lat1 < project_info.n) || (lat1 >= project_info.n && lat0 < project_info.n)) { */
- if ((lat0 >= project_info.n && lat1 <= project_info.n) || (lat1 >= project_info.n && lat0 <= project_info.n)) {
- sides[n] = 2;
- clat[n] = project_info.n;
- dlat = lat0 - lat1;
- clon[n] = (dlat == 0.0) ? lon1 : lon1 + (lon0 - lon1) * (clat[n] - lat1) / dlat;
- x_wesn_corner (&clon[n]);
- if (fabs(dlat) > 0.0 && lon_inside (clon[n], project_info.w, project_info.e)) n++;
- }
- /* if ((lon0 <= project_info.w && lon1 > project_info.w) || (lon1 <= project_info.w && lon0 > project_info.w)) { */
- if ((lon0 <= project_info.w && lon1 >= project_info.w) || (lon1 <= project_info.w && lon0 >= project_info.w)) {
- sides[n] = 3;
- clon[n] = project_info.w;
- dlon = lon0 - lon1;
- clat[n] = (dlon == 0.0) ? lat1 : lat1 + (lat0 - lat1) * (clon[n] - lon1) / dlon;
- y_wesn_corner (&clat[n]);
- if (fabs(dlon) > 0.0 && clat[n] >= project_info.s && clat[n] <= project_info.n) n++;
- }
-
- for (i = 0; i < n; i++) {
- geo_to_xy (clon[i], clat[i], &xx[i], &yy[i]);
- if (project_info.projection == POLAR && sides[i]%2) sides[i] = 4 - sides[i]; /* toggle 1 <-> 3 */
- }
-
- /* Check for corner xover */
-
- if (n < 2) return (n);
-
- if (is_wesn_corner (clon[0], clat[0])) return (1);
-
- if (is_wesn_corner (clon[1], clat[1])) {
- clon[0] = clon[1];
- clat[0] = clat[1];
- xx[0] = xx[1];
- yy[0] = yy[1];
- sides[0] = sides[1];
- return (1);
- }
-
- return (n);
- }
-
- int rect_crossing (lon0, lat0, lon1, lat1, clon, clat, xx, yy, sides)
- double lon0, lat0, lon1, lat1, *clon, *clat, *xx, *yy;
- int *sides; {
-
- /* Compute the crossover point(s) on the map boundary for rectangular projections */
- int i, n = 0;
- double x0, x1, y0, y1, d;
-
- /* Since it may not be obvious which side the line may cross, and since in some cases the two points may be
- * entirely outside the region but still cut through it, we first find all possible candidates and then decide
- * which ones are valid crossings. We may find 0, 1, or 2 intersections */
-
- geo_to_xy (lon0, lat0, &x0, &y0);
- geo_to_xy (lon1, lat1, &x1, &y1);
-
- /* First set 'almost'-corners to corners */
-
- x_rect_corner (&x0);
- x_rect_corner (&x1);
- y_rect_corner (&y0);
- y_rect_corner (&y1);
-
- /* if ((y0 > project_info.ymin && y1 <= project_info.ymin) || (y1 > project_info.ymin && y0 <= project_info.ymin)) { */
- if ((y0 >= project_info.ymin && y1 <= project_info.ymin) || (y1 >= project_info.ymin && y0 <= project_info.ymin)) {
- sides[n] = 0;
- yy[n] = project_info.ymin;
- d = y0 - y1;
- xx[n] = (d == 0.0) ? x0 : x1 + (x0 - x1) * (yy[n] - y1) / d;
- x_rect_corner (&xx[n]);
- if (xx[n] >= project_info.xmin && xx[n] <= project_info.xmax) n++;
- }
- /* if ((x0 < project_info.xmax && x1 >= project_info.xmax) || (x1 < project_info.xmax && x0 >= project_info.xmax)) { */
- if ((x0 <= project_info.xmax && x1 >= project_info.xmax) || (x1 <= project_info.xmax && x0 >= project_info.xmax)) {
- sides[n] = 1;
- xx[n] = project_info.xmax;
- d = x0 - x1;
- yy[n] = (d == 0.0) ? y0 : y1 + (y0 - y1) * (xx[n] - x1) / d;
- y_rect_corner (&yy[n]);
- if (yy[n] >= project_info.ymin && yy[n] <= project_info.ymax) n++;
- }
- /* if ((y0 < project_info.ymax && y1 >= project_info.ymax) || (y1 < project_info.ymax && y0 >= project_info.ymax)) { */
- if ((y0 <= project_info.ymax && y1 >= project_info.ymax) || (y1 <= project_info.ymax && y0 >= project_info.ymax)) {
- sides[n] = 2;
- yy[n] = project_info.ymax;
- d = y0 - y1;
- xx[n] = (d == 0.0) ? x0 : x1 + (x0 - x1) * (yy[n] - y1) / d;
- x_rect_corner (&xx[n]);
- if (xx[n] >= project_info.xmin && xx[n] <= project_info.xmax) n++;
- }
- /* if ((x0 > project_info.xmin && x1 <= project_info.xmin) || (x1 > project_info.xmin && x0 <= project_info.xmin)) { */
- if ((x0 >= project_info.xmin && x1 <= project_info.xmin) || (x1 >= project_info.xmin && x0 <= project_info.xmin)) {
- sides[n] = 3;
- xx[n] = project_info.xmin;
- d = x0 - x1;
- yy[n] = (d == 0.0) ? y0 : y1 + (y0 - y1) * (xx[n] - x1) / d;
- y_rect_corner (&yy[n]);
- if (yy[n] >= project_info.ymin && yy[n] <= project_info.ymax) n++;
- }
-
- for (i = 0; i < n; i++) xy_to_geo (&clon[i], &clat[i], xx[i], yy[i]);
-
- /* Check for corner xover */
-
- if (n < 2) return (n);
-
- if (is_rect_corner (xx[0], yy[0])) return (1);
-
- if (is_rect_corner (xx[1], yy[1])) {
- clon[0] = clon[1];
- clat[0] = clat[1];
- xx[0] = xx[1];
- yy[0] = yy[1];
- sides[0] = sides[1];
- return (1);
- }
-
- return (n);
- }
-
- int x_rect_corner (x)
- double *x; {
- if (fabs (*x) <= SMALL)
- *x = 0.0;
- else if (fabs (*x - project_info.xmax) <= SMALL)
- *x = project_info.xmax;
- }
-
- int y_rect_corner (y)
- double *y; {
- if (fabs (*y) <= SMALL)
- *y = 0.0;
- else if (fabs (*y - project_info.ymax) <= SMALL)
- *y = project_info.ymax;
- }
-
- int is_rect_corner (x, y)
- double x, y; { /* Checks if point is a corner */
- gmt_corner = -1;
- if (x == project_info.xmin) {
- if (y == project_info.ymin)
- gmt_corner = 1;
- else if (y == project_info.ymax)
- gmt_corner = 4;
- }
- else if (x == project_info.xmax) {
- if (y == project_info.ymin)
- gmt_corner = 2;
- else if (y == project_info.ymax)
- gmt_corner = 3;
- }
- return (gmt_corner > 0);
- }
-
- int x_wesn_corner (x)
- double *x; {
- /* if (fabs (fmod (fabs (*x - project_info.w), 360.0)) <= SMALL)
- *x = project_info.w;
- else if (fabs (fmod (fabs (*x - project_info.e), 360.0)) <= SMALL)
- *x = project_info.e; */
-
- if (fabs (*x - project_info.w) <= SMALL)
- *x = project_info.w;
- else if (fabs (*x - project_info.e) <= SMALL)
- *x = project_info.e;
-
- }
-
- int y_wesn_corner (y)
- double *y; {
- if (fabs (*y - project_info.s) <= SMALL)
- *y = project_info.s;
- else if (fabs (*y - project_info.n) <= SMALL)
- *y = project_info.n;
- }
-
- int is_wesn_corner (x, y)
- double x, y; { /* Checks if point is a corner */
- gmt_corner = 0;
-
- if (fmod (fabs (x - project_info.w), 360.0) == 0.0) {
- if (y == project_info.s)
- gmt_corner = 1;
- else if (y == project_info.n)
- gmt_corner = 4;
- }
- else if (fmod (fabs (x - project_info.e), 360.0) == 0.0) {
- if (y == project_info.s)
- gmt_corner = 2;
- else if (y == project_info.n)
- gmt_corner = 3;
- }
- return (gmt_corner > 0);
- }
-
- int radial_crossing (lon1, lat1, lon2, lat2, clon, clat, xx, yy, sides)
- double lon1, lat1, lon2, lat2, *clon, *clat, *xx, *yy;
- int *sides; {
- /* Computes the lon/lat of a point that is 90 spherical degrees from
- * the origin and lies on the great circle between points 1 and 2 */
-
- double dist1, dist2, delta, eps, dlon;
-
- dist1 = great_circle_dist (project_info.central_meridian, project_info.pole, lon1, lat1);
- dist2 = great_circle_dist (project_info.central_meridian, project_info.pole, lon2, lat2);
- delta = dist2 - dist1;
- eps = (delta == 0.0) ? 0.0 : (90.0 - dist1) / delta;
- dlon = lon2 - lon1;
- if (fabs (dlon) > 180.0) dlon = copysign (360.0 - fabs (dlon), -dlon);
- clon[0] = lon1 + dlon * eps;
- clat[0] = lat1 + (lat2 - lat1) * eps;
-
- geo_to_xy (clon[0], clat[0], &xx[0], &yy[0]);
- sides[0] = 1;
-
- return (1);
- }
-
- int ellipse_crossing (lon1, lat1, lon2, lat2, clon, clat, xx, yy, sides)
- double lon1, lat1, lon2, lat2, *clon, *clat, *xx, *yy;
- int *sides; {
- /* Compute the crossover point(s) on the map boundary for rectangular projections */
- int n = 0, i, jump;
- double x1, x2, y1, y2;
-
- /* Crossings here must be at the W or E borders. Lat points may only touch border */
-
- if (lat1 == -90.0) {
- sides[n] = 0;
- clon[n] = lon1;
- clat[n] = lat1;
- n = 1;
- }
- else if (lat2 == -90.0) {
- sides[n] = 0;
- clon[n] = lon2;
- clat[n] = lat2;
- n = 1;
- }
- else if (lat1 == 90.0) {
- sides[n] = 2;
- clon[n] = lon1;
- clat[n] = lat1;
- n = 1;
- }
- else if (lat2 == 90.0) {
- sides[n] = 2;
- clon[n] = lon2;
- clat[n] = lat2;
- n = 1;
- }
- else { /* May cross somewhere else */
- geo_to_xy (lon1, lat1, &x1, &y1);
- geo_to_xy (lon2, lat2, &x2, &y2);
- if ((jump = map_jump (x2, y2, x1, y1))) {
- get_crossings (xx, yy, x2, y2, x1, y1);
- if (jump == 1) { /* Add right border point first */
- d_swap (xx[0], xx[1]);
- d_swap (yy[0], yy[1]);
- }
- xy_to_geo (&clon[0], &clat[0], xx[0], yy[0]);
- xy_to_geo (&clon[1], &clat[1], xx[1], yy[1]);
- }
- n = -2; /* To signal dont change order */
- }
- if (n == 1) for (i = 0; i < n; i++) geo_to_xy (clon[i], clat[i], &xx[i], &yy[i]);
- return (n);
- }
-
- int eqdist_crossing (lon1, lat1, lon2, lat2, clon, clat, xx, yy, sides)
- double lon1, lat1, lon2, lat2, *clon, *clat, *xx, *yy;
- int *sides; {
- double angle, x, y;
-
- /* Computes the x.y of the antipole point that lies on a radius from
- * the origin through the inside point */
-
- if (eqdist_outside (lon1, lat1)) { /* Point 1 is on perimeter */
- geo_to_xy (lon2, lat2, &x, &y);
- angle = d_atan2 (y - project_info.y0, x - project_info.x0);
- xx[0] = project_info.r * cos (angle) + project_info.x0;
- yy[0] = project_info.r * sin (angle) + project_info.y0;
- clon[0] = lon1;
- clat[0] = lat1;
- }
- else { /* Point 2 is on perimeter */
- geo_to_xy (lon1, lat1, &x, &y);
- angle = d_atan2 (y - project_info.y0, x - project_info.x0);
- xx[0] = project_info.r * cos (angle) + project_info.x0;
- yy[0] = project_info.r * sin (angle) + project_info.y0;
- clon[0] = lon2;
- clat[0] = lat2;
- }
- sides[0] = 1;
-
- return (1);
- }
-
- /* Routines to add pieces of parallels or meridians */
-
- int map_path (lon1, lat1, lon2, lat2, x, y)
- double lon1, lat1, lon2, lat2;
- double *x[], *y[]; {
- if (lat1 == lat2)
- return (latpath (lat1, lon1, lon2, x, y));
- else
- return (lonpath (lon1, lat1, lat2, x, y));
- }
-
- int lonpath (lon, lat1, lat2, x, y)
- double lon, lat1, lat2;
- double *x[], *y[]; {
- int ny, n = 0, n_try, keep_trying, jump, pos;
- double dlat, *tlon, *tlat, x0, x1, y0, y1, d;
- double min_gap;
-
- min_gap = 0.1 * gmtdefs.line_step;
- if ((ny = ceil (fabs (lat2 - lat1) / gmtdefs.dlat)) == 0) return (0);
-
- ny++;
- dlat = (lat2 - lat1) / ny;
- pos = (dlat > 0.0);
-
- tlon = (double *) memory (CNULL, ny, sizeof (double), "lonpath");
- tlat = (double *) memory (CNULL, ny, sizeof (double), "lonpath");
-
- tlon[0] = lon;
- tlat[0] = lat1;
- geo_to_xy (tlon[0], tlat[0], &x0, &y0);
- while ((pos && tlat[n] < lat2) || (!pos & tlat[n] > lat2)) {
- n++;
- if (n == ny-1) {
- ny += GMT_SMALL_CHUNK;
- tlon = (double *) memory ((char *)tlon, ny, sizeof (double), "lonpath");
- tlat = (double *) memory ((char *)tlat, ny, sizeof (double), "lonpath");
- }
- n_try = 0;
- keep_trying = TRUE;
- do {
- n_try++;
- tlon[n] = lon;
- tlat[n] = tlat[n-1] + dlat;
- geo_to_xy (tlon[n], tlat[n], &x1, &y1);
- jump = map_jump (x0, y0, x1, y1) || (y0 < project_info.ymin || y0 > project_info.ymax);
- if (!jump && (d = hypot (x1 - x0, y1 - y0)) > gmtdefs.line_step)
- dlat *= 0.5;
- else if (!jump && d < min_gap)
- dlat *= 2.0;
- else
- keep_trying = FALSE;
- } while (keep_trying && n_try < 10);
- x0 = x1; y0 = y1;
- }
- tlon[n] = lon;
- tlat[n] = lat2;
- n++;
-
- if (n != ny) {
- tlon = (double *) memory ((char *)tlon, n, sizeof (double), "lonpath");
- tlat = (double *) memory ((char *)tlat, n, sizeof (double), "lonpath");
- }
-
- *x = tlon; *y = tlat;
- return (n);
- }
-
- int latpath (lat, lon1, lon2, x, y)
- double lat, lon1, lon2;
- double *x[], *y[]; {
- int nx, n = 0, n_try, keep_trying, jump, pos;
- double dlon, *tlon, *tlat, x0, x1, y0, y1, d;
- double min_gap;
-
- min_gap = 0.1 * gmtdefs.line_step;
- if ((nx = ceil (fabs (lon2 - lon1) / gmtdefs.dlon)) == 0) return (0);
-
- nx++;
- dlon = (lon2 - lon1) / nx;
- pos = (dlon > 0.0);
-
- tlon = (double *) memory (CNULL, nx, sizeof (double), "latpath");
- tlat = (double *) memory (CNULL, nx, sizeof (double), "latpath");
-
- tlon[0] = lon1;
- tlat[0] = lat;
- geo_to_xy (tlon[0], tlat[0], &x0, &y0);
- while ((pos && tlon[n] < lon2) || (!pos & tlon[n] > lon2)) {
- n++;
- if (n == nx-1) {
- nx += GMT_CHUNK;
- tlon = (double *) memory ((char *)tlon, nx, sizeof (double), "latpath");
- tlat = (double *) memory ((char *)tlat, nx, sizeof (double), "latpath");
- }
- n_try = 0;
- keep_trying = TRUE;
- do {
- n_try++;
- tlat[n] = lat;
- tlon[n] = tlon[n-1] + dlon;
- geo_to_xy (tlon[n], tlat[n], &x1, &y1);
- jump = map_jump (x0, y0, x1, y1) || (y0 < project_info.ymin || y0 > project_info.ymax);
- if (!jump && (d = hypot (x1 - x0, y1 - y0)) > gmtdefs.line_step)
- dlon *= 0.5;
- else if (!jump && d < min_gap)
- dlon *= 2.0;
- else
- keep_trying = FALSE;
- } while (keep_trying && n_try < 10);
- x0 = x1; y0 = y1;
- }
- tlon[n] = lon2;
- tlat[n] = lat;
- n++;
-
- if (n != nx) {
- tlon = (double *) memory ((char *)tlon, n, sizeof (double), "latpath");
- tlat = (double *) memory ((char *)tlat, n, sizeof (double), "latpath");
- }
-
- *x = tlon; *y = tlat;
- return (n);
- }
-
- /* Routines to do with clipping */
-
- int clip_to_map (lon, lat, np, x, y)
- double lon[], lat[];
- int np;
- double *x[], *y[]; {
- /* This routine makes sure that all points are either inside or on the map boundary
- * and returns the number of points to be used for plotting (in x,y units) */
-
- int i, n, out = 0;
- double *xx, *yy;
-
- /* First check for trivial cases: All points outside or all points inside */
-
- for (i = 0; i < np; i++) {
- map_outside (lon[i], lat[i]);
- out += (abs (gmt_x_status_new) == 2 || abs (gmt_y_status_new) == 2);
- }
- if (out == 0) { /* All points are inside map boundary */
- xx = (double *) memory (CNULL, np, sizeof (double), "clip_to_map");
- yy = (double *) memory (CNULL, np, sizeof (double), "clip_to_map");
- for (i = 0; i < np; i++) geo_to_xy (lon[i], lat[i], &xx[i], &yy[i]);
- *x = xx; *y = yy;
- n = np;
- }
- else if (out == np) /* All points are outside map boundary */
- n = 0;
- else
- n = (*map_clip) (lon, lat, np, x, y);
- return (n);
- }
-
- int rect_clip (lon, lat, n, x, y)
- double lon[], lat[];
- double *x[], *y[];
- int n; {
- int i, j = 0, k, nx, sides[2], n_alloc = GMT_CHUNK;
- double xlon[2], xlat[2], xc[2], yc[2], *xx, *yy;
-
- if (n == 0) return (0);
-
- xx = (double *) memory (CNULL, n_alloc, sizeof (double), "rect_clip");
- yy = (double *) memory (CNULL, n_alloc, sizeof (double), "rect_clip");
- map_outside (lon[0], lat[0]);
- geo_to_xy (lon[0], lat[0], &xx[0], &yy[0]);
- j += move_to_rect (xx, yy, j); /* May add 2 points, << n_alloc */
-
- /* for (i = j = 1; i < n; i++) { */
- for (i = 1; i < n; i++) {
- map_outside (lon[i], lat[i]);
- if (break_through (lon[i-1], lat[i-1], lon[i], lat[i])) {
- nx = map_crossing (lon[i-1], lat[i-1], lon[i], lat[i], xlon, xlat, xc, yc, sides);
- for (k = 0; k < nx; k++) {
- xx[j] = xc[k];
- yy[j++] = yc[k];
- if (j == (n_alloc-2)) {
- n_alloc += GMT_CHUNK;
- xx = (double *) memory ((char *)xx, n_alloc, sizeof (double), "rect_clip");
- yy = (double *) memory ((char *)yy, n_alloc, sizeof (double), "rect_clip");
- }
- }
- }
- geo_to_xy (lon[i], lat[i], &xx[j], &yy[j]);
- if (j == (n_alloc-2)) {
- n_alloc += GMT_CHUNK;
- xx = (double *) memory ((char *)xx, n_alloc, sizeof (double), "rect_clip");
- yy = (double *) memory ((char *)yy, n_alloc, sizeof (double), "rect_clip");
- }
- j += move_to_rect (xx, yy, j); /* May add 2 points, which explains the n_alloc-2 stuff */
- }
-
- xx = (double *) memory ((char *)xx, j, sizeof (double), "rect_clip");
- yy = (double *) memory ((char *)yy, j, sizeof (double), "rect_clip");
- *x = xx;
- *y = yy;
-
- return (j);
- }
-
- int wesn_clip (lon, lat, n, x, y)
- double lon[], lat[];
- double *x[], *y[];
- int n; {
- int i, j = 0, k, nx, sides[2], n_alloc = GMT_CHUNK;
- double xlon[2], xlat[2], xc[2], yc[2], *xx, *yy;
-
- if (n == 0) return (0);
-
- xx = (double *) memory (CNULL, n_alloc, sizeof (double), "wesn_clip");
- yy = (double *) memory (CNULL, n_alloc, sizeof (double), "wesn_clip");
-
- map_outside (lon[0], lat[0]);
- j = move_to_wesn (xx, yy, lon[0], lat[0], 0); /* May add 2 points, << n_alloc */
-
- for (i = 1; i < n; i++) {
- map_outside (lon[i], lat[i]);
- if (break_through (lon[i-1], lat[i-1], lon[i], lat[i])) {
- nx = map_crossing (lon[i-1], lat[i-1], lon[i], lat[i], xlon, xlat, xc, yc, sides);
- for (k = 0; k < nx; k++) {
- xx[j] = xc[k];
- yy[j++] = yc[k];
- if (j == n_alloc) {
- n_alloc += GMT_CHUNK;
- xx = (double *) memory ((char *)xx, n_alloc, sizeof (double), "wesn_clip");
- yy = (double *) memory ((char *)yy, n_alloc, sizeof (double), "wesn_clip");
- }
- }
- }
- if (j == (n_alloc-2)) {
- n_alloc += GMT_CHUNK;
- xx = (double *) memory ((char *)xx, n_alloc, sizeof (double), "wesn_clip");
- yy = (double *) memory ((char *)yy, n_alloc, sizeof (double), "wesn_clip");
- }
- j += move_to_wesn (xx, yy, lon[i], lat[i], j); /* May add 2 points, which explains the n_alloc-2 stuff */
- }
-
- xx = (double *) memory ((char *)xx, j, sizeof (double), "wesn_clip");
- yy = (double *) memory ((char *)yy, j, sizeof (double), "wesn_clip");
- *x = xx;
- *y = yy;
-
- return (j);
- }
-
- int move_to_rect (x_edge, y_edge, j)
- double x_edge[], y_edge[];
- int j; {
- int n = 0;
- double xtmp, ytmp;
-
- /* May add 0, 1, or 2 points to path */
-
- if (gmt_x_status_new == 0 && gmt_y_status_new == 0) return (1); /* Inside */
-
- if (j > 0 && gmt_x_status_new != gmt_x_status_old && gmt_y_status_new != gmt_y_status_old) { /* Must include corner */
- xtmp = x_edge[j]; ytmp = y_edge[j];
- x_edge[j] = (MIN (gmt_x_status_new, gmt_x_status_old) < 0) ? project_info.xmin : project_info.xmax;
- y_edge[j] = (MIN (gmt_y_status_new, gmt_y_status_old) < 0) ? project_info.ymin : project_info.ymax;
- j++;
- x_edge[j] = xtmp; y_edge[j] = ytmp;
- n = 1;
- }
- if (gmt_x_status_new != 0) x_edge[j] = (gmt_x_status_new < 0) ? project_info.xmin : project_info.xmax;
- if (gmt_y_status_new != 0) y_edge[j] = (gmt_y_status_new < 0) ? project_info.ymin : project_info.ymax;
- return (n + 1);
- }
-
- int move_to_wesn (x_edge, y_edge, lon, lat, j)
- double x_edge[], y_edge[], lon, lat;
- int j; {
- int n = 0;
- double xtmp, ytmp, lon_p, lat_p;
-
- /* May add 0, 1, or 2 points to path */
-
- if (j > 0 && gmt_x_status_new != gmt_x_status_old && gmt_y_status_new != gmt_y_status_old) { /* Need corner */
- xtmp = x_edge[j]; ytmp = y_edge[j];
- lon_p = (MIN (gmt_x_status_new, gmt_x_status_old) < 0) ? project_info.w : project_info.e;
- lat_p = (MIN (gmt_y_status_new, gmt_y_status_old) < 0) ? project_info.s : project_info.n;
- geo_to_xy (lon_p, lat_p, &x_edge[j], &y_edge[j]);
- j++;
- x_edge[j] = xtmp; y_edge[j] = ytmp;
- n = 1;
- }
- if (gmt_x_status_new != 0) lon = (gmt_x_status_new < 0) ? project_info.w : project_info.e;
- if (gmt_y_status_new != 0) lat = (gmt_y_status_new < 0) ? project_info.s : project_info.n;
- geo_to_xy (lon, lat, &x_edge[j], &y_edge[j]);
- return (n + 1);
- }
-
- int radial_clip (lon, lat, np, x, y)
- double lon[], lat[];
- int np;
- double *x[], *y[]; {
- int n = 0, this, i, sides[2], n_alloc = GMT_CHUNK;
- double xlon[2], xlat[2], xc[2], yc[2], xr, yr, r, scale, x0, y0, *xx, *yy;
-
- if (np == 0) return (0);
-
- xx = (double *) memory (CNULL, n_alloc, sizeof (double), "radial_clip");
- yy = (double *) memory (CNULL, n_alloc, sizeof (double), "radial_clip");
-
- if (!map_outside (lon[0], lat[0])) {
- geo_to_xy (lon[0], lat[0], &xx[0], &yy[0]);
- n++;
- }
- for (i = 1; i < np; i++) {
- this = map_outside (lon[i], lat[i]);
- if (break_through (lon[i-1], lat[i-1], lon[i], lat[i])) { /* Crossed map_boundary */
- map_crossing (lon[i-1], lat[i-1], lon[i], lat[i], xlon, xlat, xc, yc, sides);
- xx[n] = xc[0]; yy[n] = yc[0];
- n++;
- if (n == n_alloc) {
- n_alloc += GMT_CHUNK;
- xx = (double *) memory ((char *)xx, n_alloc, sizeof (double), "radial_clip");
- yy = (double *) memory ((char *)yy, n_alloc, sizeof (double), "radial_clip");
- }
- }
- geo_to_xy (lon[i], lat[i], &xr, &yr);
- if (this) { /* Project point onto perimeter */
- geo_to_xy (project_info.central_meridian, project_info.pole, &x0, &y0);
- xr -= x0; yr -= y0;
- r = hypot (xr, yr);
- scale = project_info.r / r;
- xr *= scale;
- yr *= scale;
- xr += x0; yr += y0;
- }
- xx[n] = xr; yy[n] = yr;
- n++;
- if (n == n_alloc) {
- n_alloc += GMT_CHUNK;
- xx = (double *) memory ((char *)xx, n_alloc, sizeof (double), "radial_clip");
- yy = (double *) memory ((char *)yy, n_alloc, sizeof (double), "radial_clip");
- }
- }
-
- xx = (double *) memory ((char *)xx, n, sizeof (double), "radial_clip");
- yy = (double *) memory ((char *)yy, n, sizeof (double), "radial_clip");
- *x = xx;
- *y = yy;
-
- return (n);
- }
-
- int lon_inside (lon, w, e)
- double lon, w, e; {
-
- while (lon < project_info.w) lon += 360.0;
- while (lon > project_info.e) lon -= 360.0;
-
- if (lon < w) return (FALSE);
- if (lon > e) return (FALSE);
- return (TRUE);
- }
-
- int geo_to_xy_line (lon, lat, n)
- double lon[], lat[];
- int n; {
- /* Traces the lon/lat array and returns x,y plus appropriate pen moves */
- int j, np, this, nx, sides[2], wrap = FALSE, ok;
- double xlon[2], xlat[2], xx[2], yy[2];
- double this_x, this_y, last_x, last_y, dummy[2];
-
- if (n > gmt_n_alloc) get_plot_array ();
-
- np = 0;
- geo_to_xy (lon[0], lat[0], &last_x, &last_y);
- if (!map_outside (lon[0], lat[0])) {
- gmt_x_plot[0] = last_x; gmt_y_plot[0] = last_y;
- gmt_pen[np++] = 3;
- }
- for (j = 1; j < n; j++) {
- geo_to_xy (lon[j], lat[j], &this_x, &this_y);
- this = map_outside (lon[j], lat[j]);
- nx = 0;
- if (break_through (lon[j-1], lat[j-1], lon[j], lat[j])) { /* Crossed map boundary */
- nx = map_crossing (lon[j-1], lat[j-1], lon[j], lat[j], xlon, xlat, xx, yy, sides);
- ok = ok_xovers (nx, last_x, this_x, sides, n);
- }
- if (gmt_world_map) wrap = wrap_around_check (dummy, last_x, last_y, this_x, this_y, xx, yy, sides, &nx);
- if (nx == 1) {
- gmt_x_plot[np] = xx[0]; gmt_y_plot[np] = yy[0];
- gmt_pen[np++] = pen_status ();
- if (np == gmt_n_alloc) get_plot_array ();
- }
- else if (nx == 2 && ok) {
- gmt_x_plot[np] = xx[0]; gmt_y_plot[np] = yy[0];
- gmt_pen[np++] = (wrap) ? 2 : 3;
- if (np == gmt_n_alloc) get_plot_array ();
- gmt_x_plot[np] = xx[1]; gmt_y_plot[np] = yy[1];
- gmt_pen[np++] = (wrap) ? 3 : 2;
- if (np == gmt_n_alloc) get_plot_array ();
- }
- if (!this) {
- gmt_x_plot[np] = this_x; gmt_y_plot[np] = this_y;
- gmt_pen[np++] = 2;
- if (np == gmt_n_alloc) get_plot_array ();
- }
- last_x = this_x; last_y = this_y;
- }
- return (np);
- }
-
- int ok_xovers (nx, x0, x1, sides, n)
- int nx;
- double x0, x1;
- int sides[], n; {
- if (!MAPPING) return (TRUE); /* Data is not periodic*/
- if (gmt_world_map || nx < 2) return (TRUE);
- if ((sides[0] + sides[1]) == 2) return (TRUE); /* Crossing in the n-s direction */
- if (fabs (fabs (x0 - x1) - gmt_map_width) < SMALL) return (TRUE);
- if ((sides[0] + sides[1]) != 4) return (TRUE); /* Not Crossing in the e-w direction */
- return (FALSE);
- }
-
- int compact_line (x, y, n, pen_flag, pen)
- double x[], y[];
- int n, pen[];
- BOOLEAN pen_flag; { /* TRUE if pen movements is present */
- /* compact_line will remove unnecessary points in paths */
- int i, j;
- double old_slope, new_slope, dx;
- char *flag;
-
- if (n < 3) return (n);
- flag = memory (CNULL, n, sizeof (char), "compact_line");
-
- old_slope = ( (dx = (x[1] - x[0])) == 0.0) ? copysign (MAXDOUBLE, y[1] - y[0]) : (y[1] - y[0]) / dx;
-
- for (i = 1; i < n-1; i++) {
- new_slope = ( (dx = (x[i+1] - x[i])) == 0.0) ? copysign (MAXDOUBLE, y[i+1] - y[i]) : (y[i+1] - y[i]) / dx;
- if (new_slope == old_slope && !(pen_flag && (pen[i]+pen[i+1]) > 4))
- flag[i] = 1;
- else
- old_slope = new_slope;
- }
-
- for (i = j = 1; i < n; i++) { /* i = 1 since first point must be included */
- if (flag[i] == 0) {
- x[j] = x[i];
- y[j] = y[i];
- if (pen_flag) pen[j] = pen[i];
- j++;
- }
- }
- free (flag);
-
- return (j);
- }
-
- /* Routines to transform grdfiles to/from map projections */
-
- int grdproject_init (head, x_inc, y_inc, nx, ny, dpi, offset)
- struct GRD_HEADER *head;
- double x_inc, y_inc;
- int nx, ny, dpi, offset; {
- int one;
-
- one = (offset) ? 0 : 1;
-
- if (x_inc > 0.0 && y_inc > 0.0) {
- head->nx = rint ((head->x_max - head->x_min) / x_inc) + one;
- head->ny = rint ((head->y_max - head->y_min) / y_inc) + one;
- head->x_inc = (head->x_max - head->x_min) / (head->nx - one);
- head->y_inc = (head->y_max - head->y_min) / (head->ny - one);
- }
- else if (nx > 0 && ny > 0) {
- head->nx = nx; head->ny = ny;
- head->x_inc = (head->x_max - head->x_min) / (head->nx - one);
- head->y_inc = (head->y_max - head->y_min) / (head->ny - one);
- }
- else if (dpi > 0) {
- head->nx = rint ((head->x_max - head->x_min) * dpi) + one;
- head->ny = rint ((head->y_max - head->y_min) * dpi) + one;
- head->x_inc = (head->x_max - head->x_min) / (head->nx - one);
- head->y_inc = (head->y_max - head->y_min) / (head->ny- one);
- }
- else {
- fprintf (stderr, "grdproject_init: Necessary arguments not set\n");
- exit (-1);
- }
- head->node_offset = offset;
-
- if (gmtdefs.verbose) fprintf (stderr, "grdproject: New size (nx,ny) %d by %d\n", head->nx, head->ny);
- }
-
- int grd_forward (geo, g_head, rect, r_head, max_radius, center)
- float geo[], rect[];
- struct GRD_HEADER *g_head, *r_head;
- double max_radius;
- BOOLEAN center; { /* Forward projection from geographical to rectangular grid */
- int i, j, k, ij, ii, jj, i_r, j_r, nm, di, dj, not_used = 0;
- float *weight_sum;
- double dx, dy, dr, x_0, y_0, x, y, lon, lat, delta, weight;
- double dx2, dy2, xinc2, yinc2, i_max_3r;
-
- if (project_info.projection == MERCATOR && g_head->nx == r_head->nx) {
- merc_forward (geo, g_head, rect, r_head, max_radius, center);
- return;
- }
-
- nm = r_head->nx * r_head->ny;
- weight_sum = (float *) memory (CNULL, nm, sizeof (float), "grdproject");
-
- if (max_radius == 0.0) { /* Make sensible default for search radius */
- dx = 2.0 * (r_head->x_max - r_head->x_min) / (g_head->nx);
- dy = 2.0 * (r_head->y_max - r_head->y_min) / (g_head->ny);
- if (dx < r_head->x_inc) dx = r_head->x_inc;
- if (dy < r_head->y_inc) dy = r_head->y_inc;
- max_radius = MAX (dx, dy);
- if (gmtdefs.verbose) fprintf (stderr, "grd_forward: Use max-radius = %lg\n", max_radius);
- }
- di = ceil (max_radius / r_head->x_inc);
- dj = ceil (max_radius / r_head->y_inc);
-
- dx2 = 0.5 * g_head->x_inc; dy2 = 0.5 * g_head->y_inc;
- xinc2 = 0.5 * r_head->x_inc; yinc2 = 0.5 * r_head->y_inc;
- i_max_3r = 3.0 / max_radius;
-
- for (j = ij = 0; j < g_head->ny; j++) {
- lat = g_head->y_max - j * g_head->y_inc;
- if (g_head->node_offset) lat -= dy2;
- if (project_info.projection == MERCATOR && fabs (lat) == 90.0) lat = copysign (89.99, lat);
- for (i = 0; i < g_head->nx; i++, ij++) {
- if (bad_float ((double)geo[ij])) continue;
-
- lon = g_head->x_min + i * g_head->x_inc;
- if (g_head->node_offset) lon += dx2;
- if (map_outside (lon, lat)) continue;
- geo_to_xy (lon, lat, &x_0, &y_0);
- if (center) {
- x_0 -= project_info.x0;
- y_0 -= project_info.y0;
- }
- if (r_head->node_offset) {
- ii = (x_0 == r_head->x_max) ? r_head->nx - 1 : floor ( (x_0 - r_head->x_min) / r_head->x_inc);
- jj = (y_0 == r_head->y_min) ? r_head->ny - 1 : floor ( (r_head->y_max - y_0) / r_head->y_inc);
- }
- else {
- ii = rint ( (x_0 - r_head->x_min) / r_head->x_inc);
- jj = rint ( (r_head->y_max - y_0) / r_head->y_inc);
- }
-
- for (j_r = jj - dj; j_r <= (jj + dj); j_r++) {
- if (j_r < 0 || j_r >= r_head->ny) continue;
- for (i_r = ii - di; i_r <= (ii + di); i_r++) {
- if (i_r < 0 || i_r >= r_head->nx) continue;
- k = j_r * r_head->nx + i_r;
- x = r_head->x_min + i_r * r_head->x_inc;
- y = r_head->y_max - j_r * r_head->y_inc;
- if (r_head->node_offset) {
- x += xinc2;
- y -= yinc2;
- }
- dr = hypot (x - x_0, y - y_0);
- if (dr > max_radius) continue;
- delta = dr * i_max_3r;
- weight = 1.0 / (1.0 + delta * delta);
- rect[k] += weight * geo[ij];
- weight_sum[k] += weight;
- }
- }
- }
- }
- r_head->z_min = MAXDOUBLE; r_head->z_max = -MAXDOUBLE;
- for (k = 0; k < nm; k++) {
- if (weight_sum[k] > 0.0) {
- rect[k] /= weight_sum[k];
- r_head->z_min = MIN (r_head->z_min, rect[k]);
- r_head->z_max = MAX (r_head->z_max, rect[k]);
- }
- else {
- not_used++;
- rect[k] = gmt_NaN;
- }
- }
-
- free ((char *)weight_sum);
-
- if (not_used) fprintf (stderr, "grd_forward: some projected nodes not loaded (%d)\n", not_used);
- }
-
- int grd_inverse (geo, g_head, rect, r_head, max_radius, center)
- float geo[], rect[];
- struct GRD_HEADER *g_head, *r_head;
- double max_radius;
- BOOLEAN center; { /* Transforming from rectangular projection to geographical */
- int i, j, k, ij, ii, jj, i_r, j_r, nm, di, dj, not_used = 0;
- float *weight_sum;
- double dx, dy, dr, lat_0, lon_0, x_0, y_0, lon, lat, x, y, delta, weight;
- double dx2, dy2, xinc2, yinc2, i_max_3r;
-
- if (project_info.projection == MERCATOR && g_head->nx == r_head->nx) {
- merc_inverse (geo, g_head, rect, r_head, max_radius, center);
- return;
- }
-
- nm = g_head->nx * g_head->ny;
-
- weight_sum = (float *) memory (CNULL, nm, sizeof (float), "grdproject");
-
- if (max_radius == 0.0) { /* Make sensible default for search radius */
- dx = 2.0 * (r_head->x_max - r_head->x_min) / (g_head->nx);
- dy = 2.0 * (r_head->y_max - r_head->y_min) / (g_head->ny);
- if (dx < r_head->x_inc) dx = r_head->x_inc;
- if (dy < r_head->y_inc) dy = r_head->y_inc;
- max_radius = MAX (dx, dy);
- if (gmtdefs.verbose) fprintf (stderr, "grd_inverse: Use max-radius = %lg\n", max_radius);
- }
- di = ceil (max_radius / r_head->x_inc);
- dj = ceil (max_radius / r_head->y_inc);
-
- dx2 = 0.5 * g_head->x_inc; dy2 = 0.5 * g_head->y_inc;
- xinc2 = 0.5 * r_head->x_inc; yinc2 = 0.5 * r_head->y_inc;
- i_max_3r = 3.0 / max_radius;
-
- for (j = ij = 0; j < r_head->ny; j++) {
- y_0 = r_head->y_max - j * r_head->y_inc;
- if (r_head->node_offset) y_0 -= yinc2;
- if (center) y_0 += project_info.y0;
- for (i = 0; i < r_head->nx; i++, ij++) {
- if (bad_float ((double)rect[ij])) continue;
-
- x_0 = r_head->x_min + i * r_head->x_inc;
- if (r_head->node_offset) x_0 += xinc2;
- if (center) x_0 += project_info.x0;
- xy_to_geo (&lon_0, &lat_0, x_0, y_0);
- if (g_head->node_offset) {
- ii = (lon_0 == g_head->x_max) ? g_head->nx - 1 : floor ( (lon_0 - g_head->x_min) / g_head->x_inc);
- jj = (lat_0 == g_head->y_min) ? g_head->ny - 1 : floor ( (g_head->y_max - lat_0) / g_head->y_inc);
- }
- else {
- ii = rint ( (lon_0 - g_head->x_min) / g_head->x_inc);
- jj = rint ( (g_head->y_max - lat_0) / g_head->y_inc);
- }
-
- for (j_r = jj - dj; j_r <= (jj + dj); j_r++) {
- if (j_r < 0 || j_r >= g_head->ny) continue;
- for (i_r = ii - di; i_r <= (ii + di); i_r++) {
- if (i_r < 0 || i_r >= g_head->nx) continue;
- k = j_r * g_head->nx + i_r;
- lon = g_head->x_min + i_r * g_head->x_inc;
- lat = g_head->y_max - j_r * g_head->y_inc;
- if (g_head->node_offset) {
- lon += dx2;
- lat -= dy2;
- }
- geo_to_xy (lon, lat, &x, &y);
- dr = hypot (x - x_0, y - y_0);
- if (dr > max_radius) continue;
- delta = dr * i_max_3r;
- weight = 1.0 / (1.0 + delta * delta);
- geo[k] += weight * rect[ij];
- weight_sum[k] += weight;
- }
- }
- }
- }
- g_head->z_min = MAXDOUBLE; g_head->z_max = -MAXDOUBLE;
- for (k = 0; k < nm; k++) { /* Compute weighted average */
- if (weight_sum[k] > 0.0) {
- geo[k] /= weight_sum[k];
- g_head->z_min = MIN (g_head->z_min, geo[k]);
- g_head->z_max = MAX (g_head->z_max, geo[k]);
- }
- else {
- not_used++;
- geo[k] = gmt_NaN;
- }
- }
-
- free ((char *)weight_sum);
-
- if (not_used) fprintf (stderr, "grdproject: some geographical nodes not loaded (%d)\n", not_used);
- }
-
- int merc_forward (geo, g_head, rect, r_head, max_radius, center)
- float geo[], rect[];
- struct GRD_HEADER *g_head, *r_head;
- double max_radius; /* Here max_radius is not used */
- BOOLEAN center; { /* Forward projection from geographical to mercator grid */
- int i, j, g_off, r_off;
- double dy, y, dummy, *lat_in, *lat_out, *hold, *value;
-
-
- lat_in = (double *) memory (CNULL, g_head->ny, sizeof (double), "merc_forward");
- lat_out = (double *) memory (CNULL, r_head->ny, sizeof (double), "merc_forward");
- value = (double *) memory (CNULL, r_head->ny, sizeof (double), "merc_forward");
- hold = (double *) memory (CNULL, g_head->ny, sizeof (double), "merc_forward");
-
- dy = (g_head->node_offset) ? 0.5 * g_head->y_inc : 0.0;
- for (j = 0; j < g_head->ny; j++) lat_in[j] = g_head->y_min + j * g_head->y_inc + dy;
-
- dy = (r_head->node_offset) ? 0.5 * r_head->y_inc : 0.0;
-
- for (j = 0; j < r_head->ny; j++) { /* Construct merc y-grid */
- y = r_head->y_min + j * r_head->y_inc + dy;
- if (center) y -= project_info.y0;
- xy_to_geo (&dummy, &lat_out[j], 0.0, y);
- }
-
- /* Make sure new nodes outside border are set to be on border (pixel grid only) */
-
- j = 0;
- while (j < r_head->ny && (lat_out[j] - lat_in[0]) < 0.0) lat_out[j++] = lat_in[0];
- j = r_head->ny-1;
- while (j >= 0 && (lat_out[j] - lat_in[g_head->ny-1]) > 0.0) lat_out[j--] = lat_in[g_head->ny-1];
- g_off = g_head->ny - 1;
- r_off = r_head->ny - 1;
- for (i = 0; i < r_head->nx; i++) { /* r_head->nx must == g_head->nx */
- for (j = 0; j < g_head->ny; j++) hold[g_off-j] = geo[j*g_head->nx+i]; /* Copy and reverse a column */
- intpol (lat_in, hold, g_head->ny, r_head->ny, lat_out, value, gmtdefs.interpolant);
- for (j = 0; j < r_head->ny; j++) rect[j*r_head->nx+i] = value[r_off-j]; /* Reverse and load new column */
- }
- free ((char *)lat_in);
- free ((char *)lat_out);
- free ((char *)value);
- free ((char *)hold);
- }
-
- int merc_inverse (geo, g_head, rect, r_head, max_radius, center)
- float geo[], rect[];
- struct GRD_HEADER *g_head, *r_head;
- double max_radius; /* Here max_radius is not used */
- BOOLEAN center; { /* Inverse projection from mercator to geographical grid */
- int i, j, g_off, r_off;
- double dy, y, dummy, *lat_in, *lat_out, *tmp, *val;
-
- lat_in = (double *) memory (CNULL, g_head->ny, sizeof (double), "merc_forward");
- lat_out = (double *) memory (CNULL, r_head->ny, sizeof (double), "merc_forward");
- tmp = (double *) memory (CNULL, g_head->ny, sizeof (double), "merc_forward");
- val = (double *) memory (CNULL, r_head->ny, sizeof (double), "merc_forward");
-
- dy = (g_head->node_offset) ? 0.5 * g_head->y_inc : 0.0;
- for (j = 0; j < g_head->ny; j++) lat_in[j] = g_head->y_min + j * g_head->y_inc + dy;
-
- dy = (r_head->node_offset) ? 0.5 * r_head->y_inc : 0.0;
- for (j = 0; j < r_head->ny; j++) { /* Construct merc y-grid */
- y = r_head->y_min + j * r_head->y_inc + dy;
- if (center) y -= project_info.y0;
- xy_to_geo (&dummy, &lat_out[j], 0.0, y);
- }
-
- /* Make sure new nodes outside border are set to be on border (pixel grid only) */
-
- j = 0;
- while (j < r_head->ny && (lat_in[j] - lat_out[0]) < 0.0) lat_in[j++] = lat_out[0];
- j = r_head->ny-1;
- while (j >= 0 && (lat_in[j] - lat_out[g_head->ny-1]) > 0.0) lat_in[j--] = lat_out[g_head->ny-1];
-
- g_off = g_head->ny - 1;
- r_off = r_head->ny - 1;
- for (i = 0; i < g_head->nx; i++) { /* r_head->nx must == g_head->nx */
- for (j = 0; j < r_head->ny; j++) val[r_off-j] = rect[j*r_head->nx+i]; /* Copy and reverse a column */
- intpol (lat_out, val, r_head->ny, g_head->ny, lat_in, tmp, gmtdefs.interpolant);
- for (j = 0; j < g_head->ny; j++) geo[j*g_head->nx+i] = tmp[g_off-j]; /* Copy and reverse a column */
- }
- free ((char *)lat_in);
- free ((char *)lat_out);
- free ((char *)val);
- free ((char *)tmp);
- }
-
- int two_d_to_three_d (x, y, n)
- double x[], y[];
- int n; {
- int i;
-
- /* Convert from two-D to three-D coordinates */
-
- for (i = 0; i < n; i++) xy_do_z_to_xy (x[i], y[i], project_info.z_level, &x[i], &y[i]);
- }
-
- int azim_2_angle (lon, lat, c, azim, angle)
- double lon, lat, c, azim; /* All variables in degrees */
- double *angle; {
- double lon1, lat1, x0, x1, y0, y1, sinc, cosc, sinaz, cosaz, sinl, cosl;
-
- if (project_info.projection < MERCATOR) { /* Trivial case */
- *angle = 90.0 - azim;
- return;
- }
-
- /* Find second point c spherical degrees away in the azim direction */
-
- geo_to_xy (lon, lat, &x0, &y0);
-
- azim *= D2R;
- sinaz = sin (azim);
- cosaz = cos (azim);
- c *= D2R;
- lat *= D2R;
- sinc = sin (c);
- cosc = cos (c);
- sinl = sin (lat);
- cosl = cos (lat);
-
- lon1 = lon + R2D * atan (sinc * sinaz / (cosl * cosc - sinl * sinc * cosaz));
- lat1 = R2D * d_asin (sinl * cosc + cosl * sinc * cosaz);
-
- /* Convert to x,y and get angle */
-
- geo_to_xy (lon1, lat1, &x1, &y1);
-
- *angle = d_atan2 (y1 - y0, x1 - x0) * R2D;
- }
-
- int gmt_check_R_J (clon) /* Make sure -R and -J agree for global plots; J given priority */
- double *clon; {
- double lon0;
-
- lon0 = 0.5 * (project_info.w + project_info.e);
- if (gmt_world_map && lon0 != *clon) {
- project_info.w = *clon - 180.0;
- project_info.e = *clon + 180.0;
- fprintf (stderr, "%s: GMT Warning: Central meridian set with -J (%lg) implies -R%lg/%lg/%lg/%lg\n",
- gmt_program, *clon, project_info.w, project_info.e, project_info.s, project_info.n);
- }
- else if (!gmt_world_map && !(project_info.w <= *clon && *clon <= project_info.e)) { /* Must reset meridian */
- *clon = lon0;
- fprintf (stderr, "%s: GMT Warning: Central meridian outside region, reset to %lg\n", gmt_program, lon0);
- }
- }
-
- int gmt_set_spherical () { /* Force spherical solution */
- gmtdefs.ellipsoid = N_ELLIPSOIDS - 1; /* Use equatorial radius */
- EQ_RAD = gmtdefs.ellipse[gmtdefs.ellipsoid].eq_radius;
- M_PR_DEG = PI_2 * EQ_RAD / 360.0;
- ECC = ECC2 = ECC4 = ECC6 = 0.0;
- if (gmtdefs.verbose) fprintf (stderr, "%s: GMT Warning: Spherical approximation used!\n", gmt_program);
- }
-
- int map_setinfo (xmin, xmax, ymin, ymax, scl)
- double xmin, xmax, ymin, ymax, scl; { /* Set [and rescale] parameters */
- double factor;
-
- if (project_info.gave_map_width) { /* Must rescale */
- factor = scl / ((xmax - xmin) * project_info.x_scale);
- project_info.x_scale *= factor;
- project_info.y_scale *= factor;
- }
- map_setxy (xmin, xmax, ymin, ymax);
- }
-
- int map_setxy (xmin, xmax, ymin, ymax)
- double xmin, xmax, ymin, ymax; { /* Set x/y parameters */
-
- project_info.xmax = (xmax - xmin) * project_info.x_scale;
- project_info.ymax = (ymax - ymin) * project_info.y_scale;
- project_info.x0 = -xmin * project_info.x_scale;
- project_info.y0 = -ymin * project_info.y_scale;
- }
-
- int map_clip_path (x, y, donut)
- double *x[], *y[];
- BOOLEAN *donut; {
- /* This function returns a clip path corresponding to the
- * extent of the map.
- */
-
- double *work_x, *work_y, angle, da, r0;
- int i, j, np;
-
- *donut = FALSE;
-
- if (!project_info.region) /* Rectangular map boundary */
- np = 4;
- else {
- switch (project_info.projection) {
- case LINEAR:
- case MERCATOR:
- case CYL_EQ:
- case CYL_EQDIST:
- case OBLIQUE_MERC:
- np = 4;
- break;
- case POLAR:
- *donut = (project_info.s > 0.0 && fabs (project_info.w - project_info.e) == 360.0);
- np = gmtdefs.n_lon_nodes + 1;
- if (project_info.s > 0.0) np *= 2;
- break;
- case STEREO:
- case LAMBERT:
- case LAMB_AZ_EQ:
- case ORTHO:
- case AZ_EQDIST:
- case ALBERS:
- np = 2 * gmtdefs.n_lon_nodes + 2;
- break;
- case MOLLWEIDE:
- case SINUSOIDAL:
- case ROBINSON:
- np = 2 * gmtdefs.n_lat_nodes + 2;
- break;
- case WINKEL:
- case HAMMER:
- case ECKERT:
- np = 2 * gmtdefs.n_lat_nodes + 2;
- if (project_info.s != -90.0) np += gmtdefs.n_lon_nodes - 1;
- if (project_info.n != 90.0) np += gmtdefs.n_lon_nodes - 1;
- break;
- case TM:
- case UTM:
- case CASSINI:
- np = 2 * (gmtdefs.n_lon_nodes + gmtdefs.n_lat_nodes);
- break;
- }
- }
-
- work_x = (double *)memory (CNULL, np, sizeof (double), "map_clip_on");
- work_y = (double *)memory (CNULL, np, sizeof (double), "map_clip_on");
-
- if (!project_info.region) {
- work_x[0] = work_x[3] = project_info.xmin; work_y[0] = work_y[1] = project_info.ymin;
- work_x[1] = work_x[2] = project_info.xmax; work_y[2] = work_y[3] = project_info.ymax;
-
- }
- else {
- switch (project_info.projection) { /* Fill in clip path */
- case LINEAR:
- case MERCATOR:
- case CYL_EQ:
- case CYL_EQDIST:
- case OBLIQUE_MERC:
- work_x[0] = work_x[3] = project_info.xmin; work_y[0] = work_y[1] = project_info.ymin;
- work_x[1] = work_x[2] = project_info.xmax; work_y[2] = work_y[3] = project_info.ymax;
- break;
- case LAMBERT:
- case ALBERS:
- for (i = j = 0; i <= gmtdefs.n_lon_nodes; i++, j++)
- geo_to_xy (project_info.w + i * gmtdefs.dlon, project_info.s, &work_x[j], &work_y[j]);
- for (i = 0; i <= gmtdefs.n_lon_nodes; i++, j++)
- geo_to_xy (project_info.e - i * gmtdefs.dlon, project_info.n, &work_x[j], &work_y[j]);
- break;
- case TM:
- case UTM:
- case CASSINI:
- for (i = j = 0; i < gmtdefs.n_lon_nodes; i++, j++) /* South */
- geo_to_xy (project_info.w + i * gmtdefs.dlon, project_info.s, &work_x[j], &work_y[j]);
- for (i = 0; i < gmtdefs.n_lat_nodes; j++, i++) /* East */
- geo_to_xy (project_info.e, project_info.s + i * gmtdefs.dlat, &work_x[j], &work_y[j]);
- for (i = 0; i < gmtdefs.n_lon_nodes; i++, j++) /* North */
- geo_to_xy (project_info.e - i * gmtdefs.dlon, project_info.n, &work_x[j], &work_y[j]);
- for (i = 0; i < gmtdefs.n_lat_nodes; j++, i++) /* West */
- geo_to_xy (project_info.w, project_info.n - i * gmtdefs.dlat, &work_x[j], &work_y[j]);
- break;
- case POLAR:
- r0 = project_info.r * project_info.s / project_info.n;
- if (*donut) {
- np /= 2;
- da = 2 * M_PI / np;
- for (i = 0, j = 2*np-1; i < np; i++, j--) { /* Draw outer clippath */
- angle = i * da;
- work_x[i] = project_info.r * (1.0 + cos (angle));
- work_y[i] = project_info.r * (1.0 + sin (angle));
- /* Do inner clippath and put it at end of array */
- work_x[j] = project_info.r + r0 * cos (angle);
- work_y[j] = project_info.r + r0 * sin (angle);
- }
- }
- else {
- da = fabs (project_info.e - project_info.w) / (gmtdefs.n_lon_nodes - 1);
- for (i = j = 0; i <= gmtdefs.n_lon_nodes; i++, j++) /* Draw outer clippath */
- geo_to_xy (project_info.w + i * da, project_info.n, &work_x[j], &work_y[j]);
- for (i = gmtdefs.n_lon_nodes; project_info.s > 0.0 && i >= 0; i--, j++) /* Draw inner clippath */
- geo_to_xy (project_info.w + i * da, project_info.s, &work_x[j], &work_y[j]);
- }
- break;
- case LAMB_AZ_EQ:
- case ORTHO:
- case AZ_EQDIST:
- da = 2.0 * M_PI / np;
- for (i = 0; i < np; i++) {
- angle = i * da;
- work_x[i] = project_info.r * (1.0 + cos (angle));
- work_y[i] = project_info.r * (1.0 + sin (angle));
- }
- break;
- case STEREO:
- if (project_info.polar) {
- for (i = j = 0; i <= gmtdefs.n_lon_nodes; i++, j++)
- geo_to_xy (project_info.w + i * gmtdefs.dlon, project_info.s, &work_x[j], &work_y[j]);
- for (i = 0; i <= gmtdefs.n_lon_nodes; i++, j++)
- geo_to_xy (project_info.e - i * gmtdefs.dlon, project_info.n, &work_x[j], &work_y[j]);
- }
- else {
- da = 2.0 * M_PI / np;
- for (i = 0; i < np; i++) {
- angle = i * da;
- work_x[i] = project_info.r * (1.0 + cos (angle));
- work_y[i] = project_info.r * (1.0 + sin (angle));
- }
- }
- break;
- case MOLLWEIDE:
- case SINUSOIDAL:
- case ROBINSON:
- for (i = j = 0; i <= gmtdefs.n_lat_nodes; i++, j++) /* Right */
- geo_to_xy (project_info.e, project_info.s + i * gmtdefs.dlat, &work_x[j], &work_y[j]);
- for (i = gmtdefs.n_lat_nodes; i >= 0; j++, i--) /* Left */
- geo_to_xy (project_info.w, project_info.s + i * gmtdefs.dlat, &work_x[j], &work_y[j]);
- break;
- case HAMMER:
- case WINKEL:
- case ECKERT:
- for (i = j = 0; i <= gmtdefs.n_lat_nodes; i++, j++) /* Right */
- geo_to_xy (project_info.e, project_info.s + i * gmtdefs.dlat, &work_x[j], &work_y[j]);
- for (i = 1; project_info.n != 90.0 && i < gmtdefs.n_lon_nodes; i++, j++)
- geo_to_xy (project_info.e - i * gmtdefs.dlon, project_info.n, &work_x[j], &work_y[j]);
- for (i = gmtdefs.n_lat_nodes; i >= 0; j++, i--) /* Left */
- geo_to_xy (project_info.w, project_info.s + i * gmtdefs.dlat, &work_x[j], &work_y[j]);
- for (i = 1; project_info.s != -90.0 && i < gmtdefs.n_lon_nodes; i++, j++)
- geo_to_xy (project_info.w + i * gmtdefs.dlon, project_info.s, &work_x[j], &work_y[j]);
- break;
- }
- }
-
- if (!(*donut)) np = compact_line (work_x, work_y, np, FALSE, 0);
- if (project_info.three_D) two_d_to_three_d (work_x, work_y, np);
-
- *x = work_x;
- *y = work_y;
-
- return (np);
- }
-
-