[r-cran-mapproj] 01/23: Import Upstream version 1.1-6
Andreas Tille
tille at debian.org
Fri Sep 8 08:03:05 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-mapproj.
commit bd4e4aca84fd90a0a1c10455a54e37d69e80a534
Author: Andreas Tille <tille at debian.org>
Date: Fri Sep 8 09:57:56 2017 +0200
Import Upstream version 1.1-6
---
DESCRIPTION | 9 +++
INDEX | 3 +
R/mapproj.r | 112 +++++++++++++++++++++++++
R/zzz.R | 3 +
man/map.grid.Rd | 67 +++++++++++++++
man/mapproject.Rd | 233 +++++++++++++++++++++++++++++++++++++++++++++++++++++
src/Makevars.debug | 1 +
src/aitoff.c | 31 +++++++
src/albers.c | 123 ++++++++++++++++++++++++++++
src/azequalarea.c | 24 ++++++
src/azequidist.c | 24 ++++++
src/bicentric.c | 30 +++++++
src/bonne.c | 41 ++++++++++
src/ccubrt.c | 18 +++++
src/complex.c | 89 ++++++++++++++++++++
src/conic.c | 32 ++++++++
src/cubrt.c | 35 ++++++++
src/cuts.c | 47 +++++++++++
src/cylequalarea.c | 29 +++++++
src/cylindrical.c | 24 ++++++
src/eisenlohr.c | 22 +++++
src/elco2.c | 137 +++++++++++++++++++++++++++++++
src/elliptic.c | 42 ++++++++++
src/fisheye.c | 31 +++++++
src/gall.c | 34 ++++++++
src/gilbert.c | 56 +++++++++++++
src/guyou.c | 106 ++++++++++++++++++++++++
src/harrison.c | 45 +++++++++++
src/hex.c | 129 +++++++++++++++++++++++++++++
src/homing.c | 120 +++++++++++++++++++++++++++
src/lagrange.c | 35 ++++++++
src/lambert.c | 53 ++++++++++++
src/laue.c | 29 +++++++
src/lune.c | 67 +++++++++++++++
src/map.h | 156 +++++++++++++++++++++++++++++++++++
src/mapproj_res.rc | 25 ++++++
src/mapproject.c | 158 ++++++++++++++++++++++++++++++++++++
src/mercator.c | 41 ++++++++++
src/mollweide.c | 30 +++++++
src/newyorker.c | 33 ++++++++
src/orthographic.c | 40 +++++++++
src/perspective.c | 89 ++++++++++++++++++++
src/polyconic.c | 33 ++++++++
src/rectangular.c | 27 +++++++
src/simpleconic.c | 39 +++++++++
src/sinusoidal.c | 22 +++++
src/tetra.c | 212 ++++++++++++++++++++++++++++++++++++++++++++++++
src/trapezoidal.c | 35 ++++++++
src/twocirc.c | 85 +++++++++++++++++++
src/zcoord.c | 149 ++++++++++++++++++++++++++++++++++
50 files changed, 3055 insertions(+)
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..24e2770
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,9 @@
+Package: mapproj
+Title: Map Projections
+Version: 1.1-6
+Date: 2004-03-10
+Author: Doug McIlroy. Packaged for R by Ray Brownrigg and Thomas P Minka.
+Description: Converts latitude/longitude into projected coordinates.
+Depends: maps
+License: Distribution and use for non-commercial purposes only.
+Maintainer: Thomas P Minka <surname at alum.mit.edu>
diff --git a/INDEX b/INDEX
new file mode 100644
index 0000000..e5908f4
--- /dev/null
+++ b/INDEX
@@ -0,0 +1,3 @@
+map.grid Draw a latitude/longitude grid on a projected
+ map
+mapproject Apply a Map Projection
diff --git a/R/mapproj.r b/R/mapproj.r
new file mode 100644
index 0000000..471a683
--- /dev/null
+++ b/R/mapproj.r
@@ -0,0 +1,112 @@
+"mapproject"<-
+function(x, y, projection = "", parameters = NULL, orientation = NULL)
+{
+ # minka: cleaned up handling of defaults
+ r <- NULL
+ if(!is.null(x$x)) {
+ r <- x$range[1:2]
+ y <- x$y
+ x <- x$x
+ }
+ if(length(x) != length(y))
+ stop("lengths of x and y must match")
+ if(is.null(r))
+ r <- range(x[!is.na(x)])
+ new.projection <- (projection != "")
+ if(new.projection) {
+ if(is.null(orientation)) orientation = c(90, 0, mean(r))
+ else if(length(orientation) != 3)
+ stop("orientation argument must have 3 elements")
+ }
+ else {
+ if(!exists(".Last.projection", envir = globalenv())) {
+ #stop("no previous projection")
+ return(list(x=x,y=y))
+ }
+ p <- get(".Last.projection", envir = globalenv())
+ projection <- p$projection
+ if(is.null(parameters)) parameters <- p$parameters
+ else if(length(parameters) != length(p$parameters))
+ stop(paste("expecting", length(p$parameters),
+ "parameters for", projection, "projection"))
+
+ if(is.null(orientation)) orientation <- p$orientation
+ else if(length(orientation) != 3)
+ stop("orientation argument must have 3 elements")
+ }
+ error <- .C("setproj",
+ as.character(projection),
+ as.double(parameters),
+ as.integer(length(parameters)),
+ as.double(orientation),
+ error = character(1),PACKAGE="mapproj")$error
+ if(error != "")
+ stop(error)
+ assign(".Last.projection", list(projection = projection,
+ parameters = parameters,
+ orientation = orientation),
+ envir = globalenv())
+ .C("doproj",
+ x = as.double(x),
+ y = as.double(y),
+ as.integer(length(x)),
+ range = double(4),
+ error = integer(1),
+ NAOK = TRUE,PACKAGE="mapproj")[c("x", "y", "range", "error")]
+}
+
+map.grid <- function(lim,nx=9,ny=9,labels=TRUE,pretty=TRUE,
+ cex=1,col=4,lty=2,font=2,...) {
+ # uses map.wrap from maps library
+ pretty.range <- function(lim,...) {
+ # like pretty but ensures that the range is identical:
+ # range(pretty.range(x)) == range(x)
+ x = pretty(lim,...)
+ if(abs(x[1]-lim[1]) > abs(x[2]-lim[1])) x = x[-1]
+ n = length(x)
+ if(abs(x[n]-lim[2]) > abs(x[n-1]-lim[2])) x = x[-n]
+ x[1] = lim[1]; x[length(x)] = lim[2]
+ x
+ }
+ auto.format <- function(x) {
+ # use the minimal number of digits to make x's unique
+ # similar to abbrev
+ for(digits in 0:6) {
+ s = formatC(x,digits=digits,format="f")
+ if(all(duplicated(s) == duplicated(x))) break
+ }
+ s
+ }
+ # by default, use limits of last map
+ if(missing(lim)) lim = .map.range
+ if(is.list(lim)) {
+ # first argument is a map
+ lim <- lim$range
+ }
+ if(lim[2]-lim[1] > 360) {
+ lim[2] <- lim[1] + 360
+ }
+ if(pretty) {
+ x <- pretty.range(lim[1:2],n=nx)
+ y <- pretty.range(lim[3:4],n=ny)
+ } else {
+ x <- seq(lim[1],lim[2],len=nx)
+ y <- seq(lim[3],lim[4],len=ny)
+ }
+ p = mapproject(expand.grid(x=c(seq(lim[1],lim[2],len=100),NA),y=y))
+ p = map.wrap(p)
+ lines(p,
+ col=col,lty=lty,...)
+ lines(mapproject(expand.grid(y=c(seq(lim[3],lim[4],len=100),NA),x=x)),
+ col=col,lty=lty,...)
+ if(labels) {
+ tx <- x[2]
+ xinc <- median(diff(x))
+ ty <- y[length(y)-2]
+ yinc <- median(diff(y))
+ text(mapproject(expand.grid(x=x+xinc*0.05,y=ty+yinc*0.5)),
+ labels=auto.format(x),cex=cex,adj=c(0,0),col=col,font=font,...)
+ text(mapproject(expand.grid(x=tx+xinc*0.5,y=y+yinc*0.05)),
+ labels=auto.format(y),cex=cex,adj=c(0,0),col=col,font=font,...)
+ }
+}
diff --git a/R/zzz.R b/R/zzz.R
new file mode 100644
index 0000000..36d16fb
--- /dev/null
+++ b/R/zzz.R
@@ -0,0 +1,3 @@
+.First.lib <- function(lib, pkg) {
+ library.dynam("mapproj", pkg, lib)
+}
diff --git a/man/map.grid.Rd b/man/map.grid.Rd
new file mode 100644
index 0000000..e5ddc97
--- /dev/null
+++ b/man/map.grid.Rd
@@ -0,0 +1,67 @@
+\name{map.grid}
+\alias{map.grid}
+\title{
+ Draw a latitude/longitude grid on a projected map
+}
+\description{
+ Draws a grid on an existing map.
+}
+\usage{
+map.grid(lim,nx=9,ny=9,labels=TRUE,pretty=TRUE,cex,col,lty,font,...)
+}
+\arguments{
+ \item{lim}{a vector of 4 numbers specifying
+ limits: \code{c(lon.low, lon.high, lat.low, lat.high)}.
+ \code{lim} can also be a list with a component named \code{range}, such
+ as the result of \code{map}, from which limits are taken.
+ }
+ \item{nx,ny}{the desired number of equally-spaced longitude and
+ latitude lines}
+ \item{labels}{logical to indicate if grid lines should be labeled with
+ longitude/latitude values.}
+ \item{pretty}{If \code{TRUE}, grid lines will be placed at round numbers.}
+ \item{cex}{}
+ \item{col}{}
+ \item{lty}{}
+ \item{font}{}
+ \item{...}{additional arguments passed to \code{lines} and
+ \code{text}, e.g. \code{col} to change the color of the grid and
+ \code{lty} to change the line type.}
+}
+\value{
+ Equally-spaced lines of constant longitude and lines of constant
+ latitude are superimposed on the current map, using the current projection.
+ These lines will appear curved under most projections, and give an idea of
+ how the projection works.
+}
+\seealso{\code{\link[maps]{map}}}
+\examples{
+library(maps)
+m <- map("usa",plot=FALSE)
+map("usa",project="albers",par=c(39,45))
+map.grid(m)
+
+# get unprojected world limits
+m <- map('world',plot=FALSE)
+
+# center on NYC
+map('world',proj='azequalarea',orient=c(41,-74,0))
+map.grid(m,col=2)
+points(mapproject(list(y=41,x=-74)),col=3,pch="x",cex=2)
+
+map('world',proj='orth',orient=c(41,-74,0))
+map.grid(m,col=2,nx=6,ny=5,label=FALSE,lty=2)
+points(mapproject(list(y=41,x=-74)),col=3,pch="x",cex=2)
+
+# center on Auckland
+map('world',proj='orth',orient=c(-36.92,174.6,0))
+map.grid(m,col=2,label=FALSE,lty=2)
+points(mapproject(list(y=-36.92,x=174.6)),col=3,pch="x",cex=2)
+
+m <- map('nz')
+# center on Auckland
+map('nz',proj='azequalarea',orient=c(-36.92,174.6,0))
+points(mapproject(list(y=-36.92,x=174.6)),col=3,pch="x",cex=2)
+map.grid(m,col=2)
+}
+\keyword{aplot}
diff --git a/man/mapproject.Rd b/man/mapproject.Rd
new file mode 100644
index 0000000..0694293
--- /dev/null
+++ b/man/mapproject.Rd
@@ -0,0 +1,233 @@
+\name{mapproject}
+\alias{mapproject}
+\title{
+ Apply a Map Projection
+}
+\description{
+ Converts latitude and longitude into projected coordinates.
+}
+\usage{
+mapproject(x, y, projection="", parameters=NULL, orientation=NULL)
+}
+\arguments{
+\item{x,y}{
+two vectors giving longitude and latitude coordinates
+of points on the earth's surface to be projected.
+A list containing components named \code{x} and \code{y}, giving the
+coordinates of the points to be projected may also be given.
+Missing values (\code{NA}s) are allowed.
+The coordinate system is degrees of longitude east of Greenwich
+(so the USA is bounded by negative longitudes) and degrees
+north of the equator.
+}
+\item{projection}{
+optional character string that names a map projection to
+use.
+See below for a description of this and the
+next two arguments.
+}
+\item{parameters}{
+optional numeric vector of parameters for use with the
+\code{projection} argument.
+This argument is optional only in the sense that certain
+projections do not require additional parameters.
+If a projection does require additional parameters, these
+must be given in the \code{parameters} argument.
+}
+\item{orientation}{
+ An optional vector \code{c(latitude,longitude,rotation)} which describes
+ where the "North Pole" should be when computing the projection.
+ This is mainly used for specifying the point of tangency for a planar
+ projection; you should not use it for cylindrical or conic
+ projections.
+ The third value is a clockwise rotation (in degrees).
+}}
+\value{
+list with components
+named \code{x} and \code{y}, containing the projected coordinates.
+\code{NA}s project to \code{NA}s.
+Points deemed unprojectable (such as north of 80 degrees
+latitude in the Mercator projection) are returned as \code{NA}.
+Because of the ambiguity of the first two arguments, the other
+arguments must be given by name.
+
+
+Each time \code{mapproject} is called, it leaves on frame 0 the
+dataset \code{.Last.projection}, which is a list with components \code{projection},
+\code{parameters}, and \code{orientation} giving the arguments from the
+call to \code{mapproject} or as constructed (for \code{orientation}).
+Subsequent calls to \code{mapproject} will get missing information
+from \code{.Last.projection}.
+Since \code{map} uses \code{mapproject} to do its projections, calls to
+\code{mapproject} after a call to \code{map} need not supply any arguments
+other than the data.
+}
+\details{
+Each standard projection is displayed with the Prime
+Meridian (longitude 0) being a straight vertical line, along which North
+is up.
+The orientation of nonstandard projections is specified by
+the three \code{parameters=c(lat,lon,rot)}.
+Imagine a transparent gridded sphere around the globe.
+First turn the overlay about the North Pole
+so that the Prime Meridian (longitude 0)
+of the overlay coincides with meridian \code{lon} on the globe.
+Then tilt the North Pole of the
+overlay along its Prime Meridian to latitude \code{lat} on the globe.
+Finally again turn the
+overlay about its "North Pole" so
+that its Prime Meridian coincides with the previous position
+of (the overlay) meridian \code{rot}.
+Project the desired map in
+the standard form appropriate to the overlay, but presenting
+information from the underlying globe.
+
+
+In the descriptions that follow (adapted from the McIlroy
+reference), each projection is shown as a
+function call; if it requires parameters, these are shown as
+arguments to the function.
+The descriptions are grouped into families.
+
+
+Equatorial projections centered on the Prime Meridian (longitude 0).
+Parallels are straight horizontal lines.
+\describe{
+\item{mercator()}{equally spaced straight meridians, conformal, straight
+ compass courses}
+\item{sinusoidal()}{equally spaced parallels, equal-area, same
+as \code{bonne(0)}}
+\item{cylequalarea(lat0)}{equally spaced straight meridians,
+ equal-area, true scale on \code{lat0}}
+\item{cylindrical()}{central projection on
+ tangent cylinder}
+\item{rectangular(lat0)}{equally spaced parallels, equally
+ spaced straight meridians, true scale on \code{lat0}}
+\item{gall(lat0)}{parallels
+ spaced stereographically on prime meridian, equally spaced straight
+ meridians, true scale on \code{lat0}}
+\item{mollweide()}{(homalographic) equal-area,
+ hemisphere is a circle}
+\item{gilbert()}{sphere conformally mapped on hemisphere and viewed orthographically}
+}
+
+Azimuthal projections centered on the North Pole. Parallels are
+concentric circles. Meridians are equally spaced radial lines.
+\describe{
+\item{azequidistant()}{equally spaced parallels, true distances from pole}
+\item{azequalarea()}{equal-area}
+\item{gnomonic()}{central projection on tangent plane,
+ straight great circles}
+\item{perspective(dist)}{viewed along earth's axis \code{dist}
+ earth radii from center of earth}
+\item{orthographic()}{viewed from infinity}
+\item{stereographic()}{conformal, projected from opposite pole}
+\item{laue()}{\code{radius = tan(2 * colatitude)} used in xray crystallography}
+\item{fisheye(n)}{stereographic seen through medium with refractive
+ index \code{n}}
+\item{newyorker(r)}{\code{radius = log(colatitude/r)} map from viewing
+ pedestal of radius \code{r} degrees}
+}
+
+Polar conic projections symmetric about the Prime Meridian. Parallels
+are segments of concentric circles. Except in the Bonne projection,
+meridians are equally spaced radial lines orthogonal to the parallels.
+\describe{
+\item{conic(lat0)}{central projection on cone tangent at \code{lat0}}
+\item{simpleconic(lat0,lat1)}{equally spaced parallels, true scale on \code{lat0} and \code{lat1}}
+\item{lambert(lat0,lat1)}{conformal, true scale on \code{lat0} and \code{lat1}}
+\item{albers(lat0,lat1)}{equal-area, true scale on \code{lat0} and \code{lat1}}
+\item{bonne(lat0)}{equally spaced parallels, equal-area, parallel \code{lat0}
+ developed from tangent cone}
+}
+
+Projections with bilateral symmetry about the Prime Meridian and the equator.
+\describe{
+\item{polyconic()}{parallels developed from tangent
+ cones, equally spaced along Prime Meridian}
+\item{aitoff()}{equal-area
+ projection of globe onto 2-to-1 ellipse, based on \code{azequalarea}}
+\item{lagrange()}{conformal, maps whole sphere into a circle}
+\item{bicentric(lon0)}{points plotted at true azimuth from two centers on the
+equator at longitudes \code{+lon0} and \code{-lon0}, great circles are
+straight lines (a stretched gnomonic projection)}
+\item{elliptic(lon0)}{points are
+ plotted at true distance from two centers on the equator at longitudes
+ \code{+lon0} and \code{-lon0}}
+\item{globular()}{hemisphere is circle, circular arc meridians
+ equally spaced on equator, circular arc parallels equally spaced on 0-
+ and 90-degree meridians}
+\item{vandergrinten()}{sphere is circle, meridians as
+ in \code{globular}, circular arc parallels resemble \code{mercator}}
+\item{eisenlohr()}{conformal with no singularities, shaped like polyconic}
+}
+
+Doubly periodic conformal projections.
+\describe{
+ \item{guyou}{W and E hemispheres are square}
+ \item{square}{world is square with Poles at diagonally opposite
+ corners}
+ \item{tetra}{map on tetrahedron with edge
+ tangent to Prime Meridian at S Pole,
+ unfolded into equilateral triangle}
+ \item{hex}{world is hexagon centered
+ on N Pole, N and S hemispheres are equilateral triangles}
+}
+
+Miscellaneous projections.
+\describe{
+\item{harrison(dist,angle)}{oblique
+ perspective from above the North Pole, \code{dist} earth radii from center of
+ earth, looking along the Date Line \code{angle} degrees off vertical}
+\item{trapezoidal(lat0,lat1)}{equally spaced parallels, straight meridians
+ equally spaced along parallels, true scale at \code{lat0} and \code{lat1} on Prime
+ Meridian}
+ \item{lune(lat,angle)}{conformal, polar cap above latitude \code{lat}
+ maps to convex lune with given \code{angle} at 90E and 90W}
+}
+
+Retroazimuthal projections. At every point the angle between vertical
+and a straight line to "Mecca", latitude \code{lat0} on the prime meridian,
+is the true bearing of Mecca.
+\describe{
+\item{mecca(lat0)}{equally spaced vertical meridians}
+\item{homing(lat0)}{distances to Mecca are true}
+}
+
+Maps based on the spheroid. Of geodetic quality, these projections do
+not make sense for tilted orientations.
+\describe{
+\item{sp\_mercator()}{Mercator on the spheroid.}
+\item{sp\_albers(lat0,lat1)}{Albers on the spheroid.}
+}
+} % end details
+\references{
+Richard A. Becker, and Allan R. Wilks,
+"Maps in S",
+\emph{AT\&T Bell Laboratories Statistics Research Report, 1991.}
+\url{http://www.research.att.com/areas/stat/doc/93.2.ps}
+
+M. D. McIlroy,
+documentation for
+from
+\emph{Tenth Edition UNIX Manual, Volume 1,}
+Saunders College Publishing, 1990.
+
+M. D. McIlroy, Source code for maps and map projections.
+\url{http://www.cs.dartmouth.edu/~doug/source.html}
+%\url{http://cm.bell-labs.com/cs/who/doug/map.tar.gz}
+}
+\examples{
+library(maps)
+# Bonne equal-area projection with state abbreviations
+map("state",proj='bonne', param=45)
+data(state)
+text(mapproject(state.center), state.abb)
+
+map("state",proj="albers",par=c(30,40))
+map("state",par=c(20,50)) # another Albers projection
+
+map("world",proj="gnomonic",orient=c(0,-100,0)) # example of orient
+# see map.grid for more examples
+}
+\keyword{dplot}
diff --git a/src/Makevars.debug b/src/Makevars.debug
new file mode 100644
index 0000000..97b17e7
--- /dev/null
+++ b/src/Makevars.debug
@@ -0,0 +1 @@
+PKG_CFLAGS = -Wall -pedantic -D__USE_GNU -D__USE_ISOC99
diff --git a/src/aitoff.c b/src/aitoff.c
new file mode 100644
index 0000000..467eb36
--- /dev/null
+++ b/src/aitoff.c
@@ -0,0 +1,31 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+#define Xaitwist Xaitpole.nlat
+static struct place Xaitpole;
+
+static int
+Xaitoff(struct place *place, double *x, double *y)
+{
+ struct place p;
+ copyplace(place,&p);
+ p.wlon.l /= 2.;
+ trig(&p.wlon);
+ norm(&p,&Xaitpole,&Xaitwist);
+ Xazequalarea(&p,x,y);
+ *x *= 2.;
+ return(1);
+}
+
+proj
+aitoff(void)
+{
+ latlon(0.,0.,&Xaitpole);
+ return(Xaitoff);
+}
diff --git a/src/albers.c b/src/albers.c
new file mode 100644
index 0000000..fadbddc
--- /dev/null
+++ b/src/albers.c
@@ -0,0 +1,123 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+/* For Albers formulas see Deetz and Adams "Elements of Map Projection", */
+/* USGS Special Publication No. 68, GPO 1921 */
+
+static double r0sq, r1sq, d2, n, den, sinb1, sinb2;
+static struct coord plat1, plat2;
+static southpole;
+
+static double num(double s)
+{
+ if(d2==0)
+ return(1);
+ s = d2*s*s;
+ return(1+s*(2./3+s*(3./5+s*(4./7+s*5./9))));
+}
+
+/* Albers projection for a spheroid, good only when N pole is fixed */
+
+static int
+Xspalbers(struct place *place, double *x, double *y)
+{
+ double r = sqrt(r0sq-2*(1-d2)*place->nlat.s*num(place->nlat.s)/n);
+ double t = n*place->wlon.l;
+ *y = r*cos(t);
+ *x = -r*sin(t);
+ if(!southpole)
+ *y = -*y;
+ else
+ *x = -*x;
+ return(1);
+}
+
+/* lat1, lat2: std parallels; e2: squared eccentricity */
+
+static proj albinit(double lat1, double lat2, double e2)
+{
+ double r1,r2;
+ double t;
+ for(;;) {
+ if(lat1 < -90)
+ lat1 = -180 - lat1;
+ if(lat2 > 90)
+ lat2 = 180 - lat2;
+ if(lat1 <= lat2)
+ break;
+ t = lat1; lat1 = lat2; lat2 = t;
+ }
+ if(lat2-lat1 < 1) {
+ if(lat1 > 89)
+ return(azequalarea());
+ return(0);
+ }
+ if(fabs(lat2+lat1) < 1)
+ return(cylequalarea(lat1));
+ d2 = e2;
+ den = num(1.);
+ deg2rad(lat1,&plat1);
+ deg2rad(lat2,&plat2);
+ sinb1 = plat1.s*num(plat1.s)/den;
+ sinb2 = plat2.s*num(plat2.s)/den;
+ n = (plat1.c*plat1.c/(1-e2*plat1.s*plat1.s) -
+ plat2.c*plat2.c/(1-e2*plat2.s*plat2.s)) /
+ (2*(1-e2)*den*(sinb2-sinb1));
+ r1 = plat1.c/(n*sqrt(1-e2*plat1.s*plat1.s));
+ r2 = plat2.c/(n*sqrt(2-e2*plat2.s*plat2.s));
+ r1sq = r1*r1;
+ r0sq = r1sq + 2*(1-e2)*den*sinb1/n;
+ southpole = lat1<0 && plat2.c>plat1.c;
+ return(Xspalbers);
+}
+
+proj
+sp_albers(double lat1, double lat2)
+{
+ return(albinit(lat1,lat2,EC2));
+}
+
+proj
+albers(double lat1, double lat2)
+{
+ return(albinit(lat1,lat2,0.));
+}
+
+static double scale = 1;
+static double twist = 0;
+
+void
+albscale(double x, double y, double lat, double lon)
+{
+ struct place place;
+ double alat, alon, x1,y1;
+ scale = 1;
+ twist = 0;
+ invalb(x,y,&alat,&alon);
+ twist = lon - alon;
+ deg2rad(lat,&place.nlat);
+ deg2rad(lon,&place.wlon);
+ Xspalbers(&place,&x1,&y1);
+ scale = sqrt((x1*x1+y1*y1)/(x*x+y*y));
+}
+
+void
+invalb(double x, double y, double *lat, double *lon)
+{
+ int i;
+ double sinb_den, sinp;
+ x *= scale;
+ y *= scale;
+ *lon = atan2(-x,fabs(y))/(RAD*n) + twist;
+ sinb_den = (r0sq - x*x - y*y)*n/(2*(1-d2));
+ sinp = sinb_den;
+ for(i=0; i<5; i++)
+ sinp = sinb_den/num(sinp);
+ *lat = asin(sinp)/RAD;
+}
diff --git a/src/azequalarea.c b/src/azequalarea.c
new file mode 100644
index 0000000..02a852d
--- /dev/null
+++ b/src/azequalarea.c
@@ -0,0 +1,24 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+int
+Xazequalarea(struct place *place, double *x, double *y)
+{
+ double r;
+ r = sqrt(1. - place->nlat.s);
+ *x = - r * place->wlon.s;
+ *y = - r * place->wlon.c;
+ return(1);
+}
+
+proj
+azequalarea(void)
+{
+ return(Xazequalarea);
+}
diff --git a/src/azequidist.c b/src/azequidist.c
new file mode 100644
index 0000000..ac868a3
--- /dev/null
+++ b/src/azequidist.c
@@ -0,0 +1,24 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+int
+Xazequidistant(struct place *place, double *x, double *y)
+{
+ double colat;
+ colat = PI/2 - place->nlat.l;
+ *x = -colat * place->wlon.s;
+ *y = -colat * place->wlon.c;
+ return(1);
+}
+
+proj
+azequidistant(void)
+{
+ return(Xazequidistant);
+}
diff --git a/src/bicentric.c b/src/bicentric.c
new file mode 100644
index 0000000..51a6c5c
--- /dev/null
+++ b/src/bicentric.c
@@ -0,0 +1,30 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static struct coord center;
+
+static int
+Xbicentric(struct place *place, double *x, double *y)
+{
+ if(place->wlon.c<=.01||place->nlat.c<=.01)
+ return(-1);
+ *x = -center.c*place->wlon.s/place->wlon.c;
+ *y = place->nlat.s/(place->nlat.c*place->wlon.c);
+ return(*x**x+*y**y<=9);
+}
+
+proj
+bicentric(double l)
+{
+ l = fabs(l);
+ if(l>89)
+ return(0);
+ deg2rad(l,¢er);
+ return(Xbicentric);
+}
diff --git a/src/bonne.c b/src/bonne.c
new file mode 100644
index 0000000..396f5f6
--- /dev/null
+++ b/src/bonne.c
@@ -0,0 +1,41 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static struct coord stdpar;
+static double r0;
+
+static int
+Xbonne(struct place *place, double *x, double *y)
+{
+ double r, alpha;
+ r = r0 - place->nlat.l;
+ if(r<.001)
+ if(fabs(stdpar.c)<1e-10)
+ alpha = place->wlon.l;
+ else if(fabs(place->nlat.c)==0)
+ alpha = 0;
+ else
+ alpha = place->wlon.l/(1+
+ stdpar.c*stdpar.c*stdpar.c/place->nlat.c/3);
+ else
+ alpha = place->wlon.l * place->nlat.c / r;
+ *x = - r*sin(alpha);
+ *y = - r*cos(alpha);
+ return(1);
+}
+
+proj
+bonne(double par)
+{
+ if(fabs(par*RAD) < .01)
+ return(Xsinusoidal);
+ deg2rad(par, &stdpar);
+ r0 = stdpar.c/stdpar.s + stdpar.l;
+ return(Xbonne);
+}
diff --git a/src/ccubrt.c b/src/ccubrt.c
new file mode 100644
index 0000000..abe185b
--- /dev/null
+++ b/src/ccubrt.c
@@ -0,0 +1,18 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+void
+ccubrt(double zr, double zi, double *wr, double *wi)
+{
+ double r, theta;
+ theta = atan2(zi,zr);
+ r = cubrt(hypot(zr,zi));
+ *wr = r*cos(theta/3);
+ *wi = r*sin(theta/3);
+}
diff --git a/src/complex.c b/src/complex.c
new file mode 100644
index 0000000..f3cf88b
--- /dev/null
+++ b/src/complex.c
@@ -0,0 +1,89 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+/*complex divide, defensive against overflow from
+ * * and /, but not from + and -
+ * assumes underflow yields 0.0
+ * uses identities:
+ * (a + bi)/(c + di) = ((a + bd/c) + (b - ad/c)i)/(c + dd/c)
+ * (a + bi)/(c + di) = (b - ai)/(d - ci)
+*/
+void
+cdiv(double a, double b, double c, double d, double *u, double *v)
+{
+ double r,t;
+ if(fabs(c)<fabs(d)) {
+ t = -c; c = d; d = t;
+ t = -a; a = b; b = t;
+ }
+ r = d/c;
+ t = c + r*d;
+ *u = (a + r*b)/t;
+ *v = (b - r*a)/t;
+}
+
+void
+cmul(double c1, double c2, double d1, double d2, double *e1, double *e2)
+{
+ *e1 = c1*d1 - c2*d2;
+ *e2 = c1*d2 + c2*d1;
+}
+
+void
+csq(double c1, double c2, double *e1, double *e2)
+{
+ *e1 = c1*c1 - c2*c2;
+ *e2 = c1*c2*2;
+}
+
+/* complex square root
+ * assumes underflow yields 0.0
+ * uses these identities:
+ * sqrt(x+-iy) = sqrt(r(cos(t)+-i*sin(t))
+ * = sqrt(r)(cos(t/2)+-i*sin(t/2))
+ * cos(t/2) = sin(t)/2sin(t/2) = sqrt((1+cos(t)/2)
+ * sin(t/2) = sin(t)/2cos(t/2) = sqrt((1-cos(t)/2)
+*/
+void
+map_csqrt(double c1, double c2, double *e1, double *e2)
+{
+ double r,s;
+ double x,y;
+ x = fabs(c1);
+ y = fabs(c2);
+ if(x>=y) {
+ if(x==0) {
+ *e1 = *e2 = 0;
+ return;
+ }
+ r = x;
+ s = y/x;
+ } else {
+ r = y;
+ s = x/y;
+ }
+ r *= sqrt(1+ s*s);
+ if(c1>0) {
+ *e1 = sqrt((r+c1)/2);
+ *e2 = c2/(2* *e1);
+ } else {
+ *e2 = sqrt((r-c1)/2);
+ if(c2<0)
+ *e2 = -*e2;
+ *e1 = c2/(2* *e2);
+ }
+}
+
+void map_cpow(double c1, double c2, double *d1, double *d2, double pwr)
+{
+ double theta = pwr*atan2(c2,c1);
+ double r = pow(hypot(c1,c2), pwr);
+ *d1 = r*cos(theta);
+ *d2 = r*sin(theta);
+}
diff --git a/src/conic.c b/src/conic.c
new file mode 100644
index 0000000..d6feee9
--- /dev/null
+++ b/src/conic.c
@@ -0,0 +1,32 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static struct coord stdpar;
+
+static int
+Xconic(struct place *place, double *x, double *y)
+{
+ double r;
+ if(fabs(place->nlat.l-stdpar.l) > 80.*RAD)
+ return(-1);
+ r = stdpar.c/stdpar.s - tan(place->nlat.l - stdpar.l);
+ *x = - r*sin(place->wlon.l * stdpar.s);
+ *y = - r*cos(place->wlon.l * stdpar.s);
+ if(r>3) return(0);
+ return(1);
+}
+
+proj
+conic(double par)
+{
+ if(fabs(par) <.1)
+ return(Xcylindrical);
+ deg2rad(par, &stdpar);
+ return(Xconic);
+}
diff --git a/src/cubrt.c b/src/cubrt.c
new file mode 100644
index 0000000..95543aa
--- /dev/null
+++ b/src/cubrt.c
@@ -0,0 +1,35 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+double
+cubrt(double a)
+{
+ double x,y,x1;
+ if(a==0)
+ return(0.);
+ y = 1;
+ if(a<0) {
+ y = -y;
+ a = -a;
+ }
+ while(a<1) {
+ a *= 8;
+ y /= 2;
+ }
+ while(a>1) {
+ a /= 8;
+ y *= 2;
+ }
+ x = 1;
+ do {
+ x1 = x;
+ x = (2*x1+a/(x1*x1))/3;
+ } while(fabs(x-x1)>10.e-15);
+ return(x*y);
+}
diff --git a/src/cuts.c b/src/cuts.c
new file mode 100644
index 0000000..e3b8a19
--- /dev/null
+++ b/src/cuts.c
@@ -0,0 +1,47 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+extern void abort(void);
+
+/* these routines duplicate names found in map.c. they are
+called from routines in hex.c, guyou.c, and tetra.c, which
+are in turn invoked directly from map.c. this bad organization
+arises from data hiding; only these three files know stuff
+that's necessary for the proper handling of the unusual cuts
+involved in these projections.
+
+the calling routines are not advertised as part of the library,
+and the library duplicates should never get loaded, however they
+are included to make the libary self-standing.*/
+
+int
+picut(struct place *g, struct place *og, double *cutlon)
+{
+ g; og; cutlon;
+ abort();
+ return 0;
+}
+
+int
+ckcut(struct place *g1, struct place *g2, double lon)
+{
+ g1; g2; lon;
+ abort();
+ return 0;
+}
+
+/* minka: from map.c */
+double
+reduce(double lon)
+{
+ if(lon>PI)
+ lon -= 2*PI;
+ else if(lon<-PI)
+ lon += 2*PI;
+ return(lon);
+}
diff --git a/src/cylequalarea.c b/src/cylequalarea.c
new file mode 100644
index 0000000..b1193ce
--- /dev/null
+++ b/src/cylequalarea.c
@@ -0,0 +1,29 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static double a;
+
+static int
+Xcylequalarea(struct place *place, double *x, double *y)
+{
+ *x = - place->wlon.l * a;
+ *y = place->nlat.s;
+ return(1);
+}
+
+proj
+cylequalarea(double par)
+{
+ struct coord stdp0;
+ if(par > 89.0)
+ return(0);
+ deg2rad(par, &stdp0);
+ a = stdp0.c*stdp0.c;
+ return(Xcylequalarea);
+}
diff --git a/src/cylindrical.c b/src/cylindrical.c
new file mode 100644
index 0000000..5a2fe2f
--- /dev/null
+++ b/src/cylindrical.c
@@ -0,0 +1,24 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+int
+Xcylindrical(struct place *place, double *x, double *y)
+{
+ if(fabs(place->nlat.l) > 80.*RAD)
+ return(-1);
+ *x = - place->wlon.l;
+ *y = place->nlat.s / place->nlat.c;
+ return(1);
+}
+
+proj
+cylindrical(void)
+{
+ return(Xcylindrical);
+}
diff --git a/src/eisenlohr.c b/src/eisenlohr.c
new file mode 100644
index 0000000..f2badee
--- /dev/null
+++ b/src/eisenlohr.c
@@ -0,0 +1,22 @@
+#include "map.h"
+int
+Xeisenlohr(struct place *p, double *x, double *y)
+{
+ double s1 = -sin(p->wlon.l/2);
+ double c1 = cos(p->wlon.l/2);
+ double s2 = sin(p->nlat.l/2);
+ double c2 = cos(p->nlat.l/2);
+ double t = s2/(c2+sqrt(2*p->nlat.c)*c1);
+ double c = sqrt(2/(1+t*t));
+ double q = sqrt(p->nlat.c/2);
+ double v = sqrt((c2+q*(c1+s1))/(c2+q*(c1-s1)));
+ double vi = 1/v;
+ *x = -2*log(v) + c*(v-vi);
+ *y = -2*atan(t) + c*t*(v+vi);
+ return 1;
+}
+proj
+eisenlohr()
+{
+ return Xeisenlohr;
+}
diff --git a/src/elco2.c b/src/elco2.c
new file mode 100644
index 0000000..4321073
--- /dev/null
+++ b/src/elco2.c
@@ -0,0 +1,137 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+/* elliptic integral routine, R.Bulirsch,
+ * Numerische Mathematik 7(1965) 78-90
+ * calculate integral from 0 to x+iy of
+ * (a+b*t^2)/((1+t^2)*sqrt((1+t^2)*(1+kc^2*t^2)))
+ * yields about D valid figures, where CC=10e-D
+ * for a*b>=0, except at branchpoints x=0,y=+-i,+-i/kc;
+ * there the accuracy may be reduced.
+ * fails for kc=0 or x<0
+ * return(1) for success, return(0) for fail
+ *
+ * special case a=b=1 is equivalent to
+ * standard elliptic integral of first kind
+ * from 0 to atan(x+iy) of
+ * 1/sqrt(1-k^2*(sin(t))^2) where k^2=1-kc^2
+*/
+
+#define ROOTINF 10.e18
+#define CC 1.e-6
+
+int
+elco2(double x, double y, double kc, double a, double b, double *u, double *v)
+{
+ double c,d,dn1,dn2,e,e1,e2,f,f1,f2,h,k,m,m1,m2,sy;
+ double d1[13],d2[13];
+ int i,l;
+ if(kc==0||x<0)
+ return(0);
+ sy = y>0? 1: y==0? 0: -1;
+ y = fabs(y);
+ csq(x,y,&c,&e2);
+ d = kc*kc;
+ k = 1-d;
+ e1 = 1+c;
+ cdiv2(1+d*c,d*e2,e1,e2,&f1,&f2);
+ f2 = -k*x*y*2/f2;
+ csqr(f1,f2,&dn1,&dn2);
+ if(f1<0) {
+ f1 = dn1;
+ dn1 = -dn2;
+ dn2 = -f1;
+ }
+ if(k<0) {
+ dn1 = fabs(dn1);
+ dn2 = fabs(dn2);
+ }
+ c = 1+dn1;
+ cmul(e1,e2,c,dn2,&f1,&f2);
+ cdiv(x,y,f1,f2,&d1[0],&d2[0]);
+ h = a-b;
+ d = f = m = 1;
+ kc = fabs(kc);
+ e = a;
+ a += b;
+ l = 4;
+ for(i=1;;i++) {
+ m1 = (kc+m)/2;
+ m2 = m1*m1;
+ k *= f/(m2*4);
+ b += e*kc;
+ e = a;
+ cdiv2(kc+m*dn1,m*dn2,c,dn2,&f1,&f2);
+ csqr(f1/m1,k*dn2*2/f2,&dn1,&dn2);
+ cmul(dn1,dn2,x,y,&f1,&f2);
+ x = fabs(f1);
+ y = fabs(f2);
+ a += b/m1;
+ l *= 2;
+ c = 1 +dn1;
+ d *= k/2;
+ cmul(x,y,x,y,&e1,&e2);
+ k *= k;
+
+ cmul(c,dn2,1+e1*m2,e2*m2,&f1,&f2);
+ cdiv(d*x,d*y,f1,f2,&d1[i],&d2[i]);
+ if(k<=CC)
+ break;
+ kc = sqrt(m*kc);
+ f = m2;
+ m = m1;
+ }
+ f1 = f2 = 0;
+ for(;i>=0;i--) {
+ f1 += d1[i];
+ f2 += d2[i];
+ }
+ x *= m1;
+ y *= m1;
+ cdiv2(1-y,x,1+y,-x,&e1,&e2);
+ e2 = x*2/e2;
+ d = a/(m1*l);
+ *u = atan2(e2,e1);
+ if(*u<0)
+ *u += PI;
+ a = d*sy/2;
+ *u = d*(*u) + f1*h;
+ *v = (-1-log(e1*e1+e2*e2))*a + f2*h*sy + a;
+ return(1);
+}
+
+void
+cdiv2(double c1, double c2, double d1, double d2, double *e1, double *e2)
+{
+ double t;
+ if(fabs(d2)>fabs(d1)) {
+ t = d1, d1 = d2, d2 = t;
+ t = c1, c1 = c2, c2 = t;
+ }
+ if(fabs(d1)>ROOTINF)
+ *e2 = ROOTINF*ROOTINF;
+ else
+ *e2 = d1*d1 + d2*d2;
+ t = d2/d1;
+ *e1 = (c1+t*c2)/(d1+t*d2); /* (c1*d1+c2*d2)/(d1*d1+d2*d2) */
+}
+
+/* complex square root of |x|+iy */
+void
+csqr(double c1, double c2, double *e1, double *e2)
+{
+ double r2;
+ r2 = c1*c1 + c2*c2;
+ if(r2<=0) {
+ *e1 = *e2 = 0;
+ return;
+ }
+ *e1 = sqrt((sqrt(r2) + fabs(c1))/2);
+ *e2 = c2/(*e1*2);
+}
diff --git a/src/elliptic.c b/src/elliptic.c
new file mode 100644
index 0000000..68735f8
--- /dev/null
+++ b/src/elliptic.c
@@ -0,0 +1,42 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static struct coord center;
+
+static int
+Xelliptic(struct place *place, double *x, double *y)
+{
+ double r1,r2;
+ r1 = acos(trigclamp(
+ place->nlat.c*(place->wlon.c*center.c
+ - place->wlon.s*center.s)));
+ r2 = acos(trigclamp(
+ place->nlat.c*(place->wlon.c*center.c
+ + place->wlon.s*center.s)));
+ *x = -(r1*r1 - r2*r2)/(4*center.l);
+ *y = (r1*r1+r2*r2)/2 - (center.l*center.l+*x**x);
+ if(*y < 0)
+ *y = 0;
+ *y = sqrt(*y);
+ if(place->nlat.l<0)
+ *y = -*y;
+ return(1);
+}
+
+proj
+elliptic(double l)
+{
+ l = fabs(l);
+ if(l>89)
+ return(0);
+ if(l<1)
+ return(Xazequidistant);
+ deg2rad(l,¢er);
+ return(Xelliptic);
+}
diff --git a/src/fisheye.c b/src/fisheye.c
new file mode 100644
index 0000000..c37db04
--- /dev/null
+++ b/src/fisheye.c
@@ -0,0 +1,31 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+/* refractive fisheye, not logarithmic */
+
+static double n;
+
+static int
+Xfisheye(struct place *place, double *x, double *y)
+{
+ double r;
+ double u = sin(PI/4-place->nlat.l/2)/n;
+ if(fabs(u) > .97)
+ return -1;
+ r = tan(asin(u));
+ *x = -r*place->wlon.s;
+ *y = -r*place->wlon.c;
+ return 1;
+}
+
+proj
+fisheye(double par)
+{
+ n = par;
+ return n<.1? 0: Xfisheye;
+}
diff --git a/src/gall.c b/src/gall.c
new file mode 100644
index 0000000..ec38693
--- /dev/null
+++ b/src/gall.c
@@ -0,0 +1,34 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static double scale;
+
+static int
+Xgall(struct place *place, double *x, double *y)
+{
+ /* two ways to compute tan(place->nlat.l/2) */
+ if(fabs(place->nlat.s)<.1)
+ *y = sin(place->nlat.l/2)/cos(place->nlat.l/2);
+ else
+ *y = (1-place->nlat.c)/place->nlat.s;
+ *x = -scale*place->wlon.l;
+ return 1;
+}
+
+proj
+gall(double par)
+{
+ double coshalf;
+ if(fabs(par)>80)
+ return 0;
+ par *= RAD;
+ coshalf = cos(par/2);
+ scale = cos(par)/(2*coshalf*coshalf);
+ return Xgall;
+}
diff --git a/src/gilbert.c b/src/gilbert.c
new file mode 100644
index 0000000..ca1d798
--- /dev/null
+++ b/src/gilbert.c
@@ -0,0 +1,56 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+int
+Xgilbert(struct place *p, double *x, double *y)
+{
+/* the interesting part - map the sphere onto a hemisphere */
+ struct place q;
+ q.nlat.s = tan(0.5*(p->nlat.l));
+ if(q.nlat.s > 1) q.nlat.s = 1;
+ if(q.nlat.s < -1) q.nlat.s = -1;
+ q.nlat.c = sqrt(1 - q.nlat.s*q.nlat.s);
+ q.wlon.l = p->wlon.l/2;
+ trig(&q.wlon);
+/* the dull part: present the hemisphere orthogrpahically */
+ *y = q.nlat.s;
+ *x = -q.wlon.s*q.nlat.c;
+ return(1);
+}
+
+proj
+gilbert(void)
+{
+ return(Xgilbert);
+}
+
+/* derivation of the interesting part:
+ map the sphere onto the plane by stereographic projection;
+ map the plane onto a half plane by sqrt;
+ map the half plane back to the sphere by stereographic
+ projection
+
+ n,w are original lat and lon
+ r is stereographic radius
+ primes are transformed versions
+
+ r = cos(n)/(1+sin(n))
+ r' = sqrt(r) = cos(n')/(1+sin(n'))
+
+ r'^2 = (1-sin(n')^2)/(1+sin(n')^2) = cos(n)/(1+sin(n))
+
+ this is a linear equation for sin n', with solution
+
+ sin n' = (1+sin(n)-cos(n))/(1+sin(n)+cos(n))
+
+ use standard formula: tan x/2 = (1-cos x)/sin x = sin x/(1+cos x)
+ to show that the right side of the last equation is tan(n/2)
+*/
+
+
diff --git a/src/guyou.c b/src/guyou.c
new file mode 100644
index 0000000..fbd25cb
--- /dev/null
+++ b/src/guyou.c
@@ -0,0 +1,106 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static struct place gywhem, gyehem;
+static struct coord gytwist;
+static double gyconst, gykc, gyside;
+
+
+static void
+dosquare(double z1, double z2, double *x, double *y)
+{
+ double w1,w2;
+ w1 = z1 -1;
+ if(fabs(w1*w1+z2*z2)>.000001) {
+ cdiv(z1+1,z2,w1,z2,&w1,&w2);
+ w1 *= gyconst;
+ w2 *= gyconst;
+ if(w1<0)
+ w1 = 0;
+ elco2(w1,w2,gykc,1.,1.,x,y);
+ } else {
+ *x = gyside;
+ *y = 0;
+ }
+}
+
+int
+Xguyou(struct place *place, double *x, double *y)
+{
+ int ew; /*which hemisphere*/
+ double z1,z2;
+ struct place pl;
+ ew = place->wlon.l<0;
+ copyplace(place,&pl);
+ norm(&pl,ew?&gyehem:&gywhem,&gytwist);
+ Xstereographic(&pl,&z1,&z2);
+ dosquare(z1/2,z2/2,x,y);
+ if(!ew)
+ *x -= gyside;
+ return(1);
+}
+
+proj
+guyou(void)
+{
+ double junk;
+ gykc = 1/(3+2*sqrt(2.));
+ gyconst = -(1+sqrt(2.));
+ elco2(-gyconst,0.,gykc,1.,1.,&gyside,&junk);
+ gyside *= 2;
+ latlon(0.,90.,&gywhem);
+ latlon(0.,-90.,&gyehem);
+ deg2rad(0.,&gytwist);
+ return(Xguyou);
+}
+
+int
+guycut(struct place *g, struct place *og, double *cutlon)
+{
+ int c;
+ c = picut(g,og,cutlon);
+ if(c!=1)
+ return(c);
+ *cutlon = 0.;
+ if(g->nlat.c<.7071||og->nlat.c<.7071)
+ return(ckcut(g,og,0.));
+ return(1);
+}
+
+static int
+Xsquare(struct place *place, double *x, double *y)
+{
+ double z1,z2;
+ double r, theta;
+ struct place p;
+ copyplace(place,&p);
+ if(place->nlat.l<0) {
+ p.nlat.l = -p.nlat.l;
+ p.nlat.s = -p.nlat.s;
+ }
+ if(p.nlat.l<FUZZ && fabs(p.wlon.l)>PI-FUZZ){
+ *y = -gyside/2;
+ *x = p.wlon.l>0?0:gyside;
+ return(1);
+ }
+ Xstereographic(&p,&z1,&z2);
+ r = sqrt(sqrt(hypot(z1,z2)/2));
+ theta = atan2(z1,-z2)/4;
+ dosquare(r*sin(theta),-r*cos(theta),x,y);
+ if(place->nlat.l<0)
+ *y = -gyside - *y;
+ return(1);
+}
+
+proj
+square(void)
+{
+ guyou();
+ return(Xsquare);
+}
diff --git a/src/harrison.c b/src/harrison.c
new file mode 100644
index 0000000..c0494ab
--- /dev/null
+++ b/src/harrison.c
@@ -0,0 +1,45 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static double v3,u2,u3,a,b; /*v=view,p=obj,u=unit.y*/
+
+static int
+Xharrison(struct place *place, double *x, double *y)
+{
+ double p1 = -place->nlat.c*place->wlon.s;
+ double p2 = -place->nlat.c*place->wlon.c;
+ double p3 = place->nlat.s;
+ double d = b + u3*p2 - u2*p3;
+ double t;
+ if(d < .01)
+ return -1;
+ t = a/d;
+ if(v3*place->nlat.s < 1.)
+ return -1;
+ *y = t*p2*u2 + (v3-t*(v3-p3))*u3;
+ *x = t*p1;
+ if(t < 0)
+ return 0;
+ if(*x * *x + *y * *y > 16)
+ return -1;
+ return 1;
+}
+
+proj
+harrison(double r, double alpha)
+{
+ u2 = cos(alpha*RAD);
+ u3 = sin(alpha*RAD);
+ v3 = r;
+ b = r*u2;
+ a = 1 + b;
+ if(r<1.001 || a<sqrt(r*r-1))
+ return 0;
+ return Xharrison;
+}
diff --git a/src/hex.c b/src/hex.c
new file mode 100644
index 0000000..6f5bebc
--- /dev/null
+++ b/src/hex.c
@@ -0,0 +1,129 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+#define BIG 1.e15
+#define HFUZZ .0001
+
+static double hcut[3] ;
+static double kr[3] = { .5, -1., .5 };
+static double ki0[3] = { -1., 0., 1. };
+static double ki[3]; /*ki0*sqrt(3)/2*/
+static double cr[3];
+static double ci[3];
+static struct place hem;
+static struct coord twist;
+static double rootroot3, hkc;
+static double w2;
+static double rootk;
+
+static void
+reflect(int i, double wr, double wi, double *x, double *y)
+{
+ double pr,pi,l;
+ pr = cr[i]-wr;
+ pi = ci[i]-wi;
+ l = 2*(kr[i]*pr + ki[i]*pi);
+ *x = wr + l*kr[i];
+ *y = wi + l*ki[i];
+}
+
+static int
+Xhex(struct place *place, double *x, double *y)
+{
+ int ns;
+ int i;
+ double zr,zi;
+ double sr,si,tr,ti,ur,ui,vr,vi,yr,yi;
+ struct place p;
+ copyplace(place,&p);
+ ns = (place->nlat.l >= 0);
+ if(!ns) {
+ p.nlat.l = -p.nlat.l;
+ p.nlat.s = -p.nlat.s;
+ }
+ if(p.nlat.l<HFUZZ) {
+ for(i=0;i<3;i++) {
+ if(fabs(reduce(p.wlon.l-hcut[i]))<HFUZZ) {
+ if(i==2) {
+ *x = 2*cr[0] - cr[1];
+ *y = 0;
+ } else {
+ *x = cr[1];
+ *y = 2*ci[2*i];
+ }
+ return(1);
+ }
+ }
+ p.nlat.l = HFUZZ;
+ trig(&p.nlat);
+ }
+ norm(&p,&hem,&twist);
+ Xstereographic(&p,&zr,&zi);
+ zr /= 2;
+ zi /= 2;
+ cdiv(1-zr,-zi,1+zr,zi,&sr,&si);
+ csq(sr,si,&tr,&ti);
+ ccubrt(1+3*tr,3*ti,&ur,&ui);
+ map_csqrt(ur-1,ui,&vr,&vi);
+ cdiv(rootroot3+vr,vi,rootroot3-vr,-vi,&yr,&yi);
+ yr /= rootk;
+ yi /= rootk;
+ elco2(fabs(yr),yi,hkc,1.,1.,x,y);
+ if(yr < 0)
+ *x = w2 - *x;
+ if(!ns) reflect(hcut[0]>place->wlon.l?0:
+ hcut[1]>=place->wlon.l?1:
+ 2,*x,*y,x,y);
+ return(1);
+}
+
+proj
+map_hex(void)
+{
+ int i;
+ double t;
+ double root3;
+ double c,d;
+ struct place p;
+ hcut[2] = PI;
+ hcut[1] = hcut[2]/3;
+ hcut[0] = -hcut[1];
+ root3 = sqrt(3.);
+ rootroot3 = sqrt(root3);
+ t = 15 -8*root3;
+ hkc = t*(1-sqrt(1-1/(t*t)));
+ elco2(BIG,0.,hkc,1.,1.,&w2,&t);
+ w2 *= 2;
+ rootk = sqrt(hkc);
+ latlon(90.,90.,&hem);
+ latlon(90.,0.,&p);
+ Xhex(&p,&c,&t);
+ latlon(0.,0.,&p);
+ Xhex(&p,&d,&t);
+ for(i=0;i<3;i++) {
+ ki[i] = ki0[i]*root3/2;
+ cr[i] = c + (c-d)*kr[i];
+ ci[i] = (c-d)*ki[i];
+ }
+ deg2rad(0.,&twist);
+ return(Xhex);
+}
+
+int
+hexcut(struct place *g, struct place *og, double *cutlon)
+{
+ int t,i;
+ if(g->nlat.l>=-HFUZZ&&og->nlat.l>=-HFUZZ)
+ return(1);
+ for(i=0;i<3;i++) {
+ t = ckcut(g,og,*cutlon=hcut[i]);
+ if(t!=1) return(t);
+ }
+ return(1);
+}
diff --git a/src/homing.c b/src/homing.c
new file mode 100644
index 0000000..5912134
--- /dev/null
+++ b/src/homing.c
@@ -0,0 +1,120 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static struct coord p0; /* standard parallel */
+
+int first;
+
+static struct coord az; /* azimuth of p0 seen from place */
+static struct coord rad; /* angular dist from place to p0 */
+
+static int
+azimuth(struct place *place)
+{
+ if(place->nlat.c < FUZZ) {
+ az.l = PI/2 + place->nlat.l - place->wlon.l;
+ trig(&az);
+ rad.l = fabs(place->nlat.l - p0.l);
+ if(rad.l > PI)
+ rad.l = 2*PI - rad.l;
+ trig(&rad);
+ return 1;
+ }
+ rad.c = trigclamp(p0.s*place->nlat.s + /* law of cosines */
+ p0.c*place->nlat.c*place->wlon.c);
+ rad.s = sqrt(1 - rad.c*rad.c);
+ if(fabs(rad.s) < .001) {
+ az.s = 0;
+ az.c = 1;
+ } else {
+ az.s = trigclamp(p0.c*place->wlon.s/rad.s); /* sines */
+ az.c = trigclamp((p0.s - rad.c*place->nlat.s)
+ /(rad.s*place->nlat.c));
+ }
+ rad.l = atan2(rad.s, rad.c);
+ return 1;
+}
+
+static int
+Xmecca(struct place *place, double *x, double *y)
+{
+ if(!azimuth(place))
+ return 0;
+ *x = -place->wlon.l;
+ *y = fabs(az.s)<.02? -az.c*rad.s/p0.c: *x*az.c/az.s;
+ return fabs(*y)>2? -1:
+ rad.c<0? 0:
+ 1;
+}
+
+proj
+mecca(double par)
+{
+ first = 1;
+ if(fabs(par)>80.)
+ return(0);
+ deg2rad(par,&p0);
+ return(Xmecca);
+}
+
+static int
+Xhoming(struct place *place, double *x, double *y)
+{
+ if(!azimuth(place))
+ return 0;
+ *x = -rad.l*az.s;
+ *y = -rad.l*az.c;
+ return place->wlon.c<0? 0: 1;
+}
+
+proj
+homing(double par)
+{
+ first = 1;
+ if(fabs(par)>80.)
+ return(0);
+ deg2rad(par,&p0);
+ return(Xhoming);
+}
+
+int
+hlimb(double *lat, double *lon, double res)
+{
+ if(first) {
+ *lon = -90;
+ *lat = -90;
+ first = 0;
+ return 0;
+ }
+ *lat += res;
+ if(*lat <= 90)
+ return 1;
+ if(*lon == 90)
+ return -1;
+ *lon = 90;
+ *lat = -90;
+ return 0;
+}
+
+int
+mlimb(double *lat, double *lon, double res)
+{
+ int ret = !first;
+ if(fabs(p0.s) < .01)
+ return -1;
+ if(first) {
+ *lon = -180;
+ first = 0;
+ } else
+ *lon += res;
+ if(*lon > 180)
+ return -1;
+ *lat = atan(-cos(*lon*RAD)/p0.s*p0.c)/RAD;
+ return ret;
+}
diff --git a/src/lagrange.c b/src/lagrange.c
new file mode 100644
index 0000000..2fde4fe
--- /dev/null
+++ b/src/lagrange.c
@@ -0,0 +1,35 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static int
+Xlagrange(struct place *place, double *x, double *y)
+{
+ double z1,z2;
+ double w1,w2,t1,t2;
+ struct place p;
+ copyplace(place,&p);
+ if(place->nlat.l<0) {
+ p.nlat.l = -p.nlat.l;
+ p.nlat.s = -p.nlat.s;
+ }
+ Xstereographic(&p,&z1,&z2);
+ map_csqrt(-z2/2,z1/2,&w1,&w2);
+ cdiv(w1-1,w2,w1+1,w2,&t1,&t2);
+ *y = -t1;
+ *x = t2;
+ if(place->nlat.l<0)
+ *y = -*y;
+ return(1);
+}
+
+proj
+lagrange(void)
+{
+ return(Xlagrange);
+}
diff --git a/src/lambert.c b/src/lambert.c
new file mode 100644
index 0000000..7fb3e0d
--- /dev/null
+++ b/src/lambert.c
@@ -0,0 +1,53 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static struct coord stdp0, stdp1;
+static double k;
+
+static int
+Xlambert(struct place *place, double *x, double *y)
+{
+ double r;
+ if(place->nlat.l < -80.*RAD)
+ return(-1);
+ if(place->nlat.l > 89.*RAD)
+ r = 0; /* slovenly */
+ else
+ r = stdp0.c*exp(0.5*k*log(
+ (1+stdp0.s)*(1-place->nlat.s)/((1-stdp0.s)*(1+place->nlat.s))));
+ if(stdp1.l<0.)
+ r = -r;
+ *x = - r*sin(k * place->wlon.l);
+ *y = - r*cos(k * place->wlon.l);
+ return(1);
+}
+
+proj
+lambert(double par0, double par1)
+{
+ double temp;
+ if(fabs(par0)>fabs(par1)){
+ temp = par0;
+ par0 = par1;
+ par1 = temp;
+ }
+ deg2rad(par0, &stdp0);
+ deg2rad(par1, &stdp1);
+ if(fabs(par1+par0)<.1)
+ return(mercator());
+ if(fabs(par0)>89.5||fabs(par1)>89.5)
+ return(0);
+ if(fabs(par1-par0)<.1)
+ /* series expansion about stdp1.s = stdp0.s */
+ k = stdp0.s + 0.5*(stdp1.s - stdp0.s);
+ else
+ k = 2*log(stdp1.c/stdp0.c)/log(
+ (1+stdp0.s)*(1-stdp1.s)/((1-stdp0.s)*(1+stdp1.s)));
+ return(Xlambert);
+}
diff --git a/src/laue.c b/src/laue.c
new file mode 100644
index 0000000..e17ae0f
--- /dev/null
+++ b/src/laue.c
@@ -0,0 +1,29 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+
+static int
+Xlaue(struct place *place, double *x, double *y)
+{
+ double r;
+ if(place->nlat.l<PI/4+FUZZ)
+ return(-1);
+ r = tan(PI-2*place->nlat.l);
+ if(r>3)
+ return(-1);
+ *x = - r * place->wlon.s;
+ *y = - r * place->wlon.c;
+ return(1);
+}
+
+proj
+laue(void)
+{
+ return(Xlaue);
+}
diff --git a/src/lune.c b/src/lune.c
new file mode 100644
index 0000000..6383a0b
--- /dev/null
+++ b/src/lune.c
@@ -0,0 +1,67 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+int Xstereographic(struct place *place, double *x, double *y);
+
+static struct place eastpole;
+static struct place westpole;
+static double eastx, easty;
+static double westx, westy;
+static double scale;
+static double pwr;
+
+/* conformal map w = ((1+z)^A - (1-z)^A)/((1+z)^A + (1-z)^A),
+ where A<1, maps unit circle onto a convex lune with x= +-1
+ mapping to vertices of angle A*PI at w = +-1 */
+
+/* there are cuts from E and W poles to S pole,
+ in absence of a cut routine, error is returned for
+ points outside a polar cap through E and W poles */
+
+static int Xlune(struct place *place, double *x, double *y)
+{
+ double stereox, stereoy;
+ double z1x, z1y, z2x, z2y;
+ double w1x, w1y, w2x, w2y;
+ double numx, numy, denx, deny;
+ if(place->nlat.l < eastpole.nlat.l-FUZZ)
+ return -1;
+ Xstereographic(place, &stereox, &stereoy);
+ stereox *= scale;
+ stereoy *= scale;
+ z1x = 1 + stereox;
+ z1y = stereoy;
+ z2x = 1 - stereox;
+ z2y = -stereoy;
+ map_cpow(z1x,z1y,&w1x,&w1y,pwr);
+ map_cpow(z2x,z2y,&w2x,&w2y,pwr);
+ numx = w1x - w2x;
+ numy = w1y - w2y;
+ denx = w1x + w2x;
+ deny = w1y + w2y;
+ cdiv(numx, numy, denx, deny, x, y);
+ return 1;
+}
+
+proj
+lune(double lat, double theta)
+{
+ deg2rad(lat, &eastpole.nlat);
+ deg2rad(-90.,&eastpole.wlon);
+ deg2rad(lat, &westpole.nlat);
+ deg2rad(90. ,&westpole.wlon);
+ Xstereographic(&eastpole, &eastx, &easty);
+ Xstereographic(&westpole, &westx, &westy);
+ if(fabs(easty)>FUZZ || fabs(westy)>FUZZ ||
+ fabs(eastx+westx)>FUZZ)
+ abort();
+ scale = 1/eastx;
+ pwr = theta/180;
+ return Xlune;
+}
diff --git a/src/map.h b/src/map.h
new file mode 100644
index 0000000..da32104
--- /dev/null
+++ b/src/map.h
@@ -0,0 +1,156 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+
+#ifndef PI
+#define PI 3.1415926535897932384626433832795028841971693993751
+#endif
+
+#define TWOPI (2*PI)
+#define RAD (PI/180)
+double hypot(double, double); /* sqrt(a*a+b*b) */
+double tan(double); /* not in K&R library */
+
+#define ECC .08227185422 /* eccentricity of earth */
+#define EC2 .006768657997
+
+#define FUZZ .0001
+#define UNUSED 0.0 /* a dummy double parameter */
+
+struct coord {
+ double l; /* lat or lon in radians*/
+ double s; /* sin */
+ double c; /* cos */
+};
+struct place {
+ struct coord nlat;
+ struct coord wlon;
+};
+
+typedef int (*proj)(struct place *, double *, double *);
+
+struct pindex { /* index of known projections */
+ char *name; /* name of projection */
+ proj (*prog)(double, double);
+ /* pointer to projection function */
+ int npar; /* number of params */
+ int (*cut)(struct place *, struct place *, double *);
+ /* function that handles cuts--eg longitude 180 */
+ int poles; /*1 S pole is a line, 2 N pole is, 3 both*/
+ int spheroid; /* poles must be at 90 deg if nonzero */
+ int (*limb)(double *lat, double *lon, double resolution);
+ /* get next place on limb */
+ /* return -1 if done, 0 at gap, else 1 */
+};
+
+
+proj aitoff(void);
+proj albers(double, double);
+int Xazequalarea(struct place *, double *, double *);
+proj azequalarea(void);
+int Xazequidistant(struct place *, double *, double *);
+proj azequidistant(void);
+proj bicentric(double);
+proj bonne(double);
+proj conic(double);
+proj cylequalarea(double);
+int Xcylindrical(struct place *, double *, double *);
+proj cylindrical(void);
+proj elliptic(double);
+proj eisenlohr(void);
+proj fisheye(double);
+proj gall(double);
+proj gilbert(void);
+proj globular(void);
+proj gnomonic(void);
+int guycut(struct place *, struct place *, double *);
+int Xguyou(struct place *, double *, double *);
+proj guyou(void);
+proj harrison(double, double);
+proj map_hex(void);
+proj homing(double);
+int hlimb(double*, double*, double resolution);
+proj lagrange(void);
+proj lambert(double, double);
+proj laue(void);
+proj lune(double, double);
+proj loxodromic(double); /* not in library */
+proj mecca(double);
+int mlimb(double*, double*, double resolution);
+proj mercator(void);
+proj mollweide(void);
+proj newyorker(double);
+proj ortelius(double, double); /* not in library */
+int Xorthographic(struct place *place, double *x, double *y);
+proj orthographic(void);
+int olimb(double*, double*, double);
+proj map_perspective(double);
+int plimb(double*, double*, double resolution);
+int Xpolyconic(struct place *, double *, double *);
+proj polyconic(void);
+proj rectangular(double);
+proj simpleconic(double, double);
+int Xsinusoidal(struct place *, double *, double *);
+proj sinusoidal(void);
+proj sp_albers(double, double);
+proj sp_mercator(void);
+proj square(void);
+int Xstereographic(struct place *, double *, double *);
+proj stereographic(void);
+int Xtetra(struct place *, double *, double *);
+int tetracut(struct place *, struct place *, double *);
+proj tetra(void);
+proj trapezoidal(double, double);
+proj vandergrinten(void);
+proj wreath(double, double); /* not in library */
+
+void findxy(double, double *, double *);
+void albscale(double, double, double, double);
+void invalb(double, double, double *, double *);
+
+void cdiv(double, double, double, double, double *, double *);
+void cmul(double, double, double, double, double *, double *);
+void map_cpow(double, double, double *, double *, double);
+void csq(double, double, double *, double *);
+void map_csqrt(double, double, double *, double *);
+void ccubrt(double, double, double *, double *);
+double cubrt(double);
+int elco2(double, double, double, double, double, double *, double *);
+void cdiv2(double, double, double, double, double *, double *);
+void csqr(double, double, double *, double *);
+
+void orient(double, double, double);
+void latlon(double, double, struct place *);
+void deg2rad(double, struct coord *);
+double trigclamp(double);
+void trig(struct coord *);
+void normalize(struct place *);
+void invert(struct place *);
+void norm(struct place *, struct place *, struct coord *);
+void printp(struct place *);
+void copyplace(struct place *, struct place *);
+
+int picut(struct place *, struct place *, double *);
+int ckcut(struct place *, struct place *, double);
+double reduce(double);
+
+/* (minka)
+void getsyms(char *);
+int putsym(struct place *, char *, double, int);
+void filerror(char *s, char *f);
+void error(char *s);
+int doproj(struct place *, int *, int *);
+int cpoint(int, int, int);
+int plotpt(struct place *, int);
+int nocut(struct place *, struct place *, double *);
+*/
+
+extern int (*projection)(struct place *, double *, double *);
diff --git a/src/mapproj_res.rc b/src/mapproj_res.rc
new file mode 100644
index 0000000..d47a8a0
--- /dev/null
+++ b/src/mapproj_res.rc
@@ -0,0 +1,25 @@
+#include <windows.h>
+#include "Rversion.h"
+
+VS_VERSION_INFO VERSIONINFO
+FILEVERSION R_FILEVERSION
+PRODUCTVERSION 3,0,0,0
+FILEFLAGSMASK 0x3L
+FILEOS VOS__WINDOWS32
+FILETYPE VFT_APP
+BEGIN
+ BLOCK "StringFileInfo"
+ BEGIN
+ BLOCK "040904E4"
+ BEGIN
+ VALUE "FileDescription", "DLL for R package `mapproj'\0"
+ VALUE "FileVersion", "1.1\0"
+ VALUE "Compiled under R Version", R_MAJOR "." R_MINOR " (" R_YEAR "-" R_MONTH "-" R_DAY ")\0"
+ VALUE "Project info", "http://www.r-project.org\0"
+ END
+ END
+ BLOCK "VarFileInfo"
+ BEGIN
+ VALUE "Translation", 0x409, 1252
+ END
+END
diff --git a/src/mapproject.c b/src/mapproject.c
new file mode 100644
index 0000000..86e8cf9
--- /dev/null
+++ b/src/mapproject.c
@@ -0,0 +1,158 @@
+#include "R.h"
+#include "map.h"
+
+#define MIN(a,b) (a)<(b)?(a):(b)
+#define MAX(a,b) (a)>(b)?(a):(b)
+#define ABS(x) ((x)<0?-(x):(x))
+#define RAD2DEG(x) ((x)*180/PI)
+#define DEG2RAD(x) ((x)*PI/180)
+
+/* must be consistent with map.h in maps library */
+#define XMIN 0
+#define XMAX 1
+#define YMIN 2
+#define YMAX 3
+
+static int (*projfun)();
+
+struct index {
+ char *name;
+ proj (*prog)();
+ int npar;
+} mapindex[] = {
+ {"mercator", mercator, 0},
+ {"cylindrical", cylindrical, 0},
+ {"cylequalarea", cylequalarea, 1},
+ {"rectangular", rectangular, 1},
+ {"trapezoidal", trapezoidal, 2},
+ {"lune",lune,2},
+ {"gall", gall, 1},
+ {"sinusoidal", sinusoidal, 0},
+ {"mollweide", mollweide, 0},
+ {"gilbert", gilbert, 0},
+ {"azequidistant", azequidistant, 0},
+ {"azequalarea", azequalarea, 0},
+ {"gnomonic", gnomonic, 0},
+ {"perspective", map_perspective, 1},
+ {"harrison", harrison, 2},
+ {"orthographic", orthographic, 0},
+ {"stereographic", stereographic, 0},
+ {"laue", laue, 0},
+ {"fisheye", fisheye, 1},
+ {"newyorker", newyorker, 1},
+ {"conic", conic, 1},
+ {"lambert", lambert, 2},
+ {"albers", albers, 2},
+ {"bonne", bonne, 1},
+ {"polyconic", polyconic, 0},
+ {"aitoff", aitoff, 0},
+ {"globular", globular, 0},
+ {"vandergrinten", vandergrinten, 0},
+ {"eisenlohr", eisenlohr, 0},
+ {"guyou",guyou,0},
+ {"square",square,0},
+ {"tetra",tetra,0},
+ {"hex",map_hex,0},
+ {"lagrange",lagrange,0},
+ {"bicentric", bicentric, 1},
+ {"elliptic", elliptic, 1},
+ {"mecca", mecca, 1},
+ {"simpleconic", simpleconic, 2},
+ {"homing", homing, 1},
+ {"sp_mercator", sp_mercator, 0},
+ {"sp_albers", sp_albers, 2},
+ {NULL,NULL,0},
+};
+
+void setproj(name, par, n, o, error)
+ char **name, **error;
+ double par[], o[];
+ long *n;
+{
+ struct index *i, *theproj = 0;
+ static char errbuf[200];
+
+ *error = "";
+ if(**name == 0) {
+ *error = "Null projection specified";
+ return;
+ }
+ for(i = mapindex; i->name != 0; i++) {
+ if(strncmp(*name, i->name, strlen(*name)) == 0) {
+ if(theproj) {
+ sprintf(errbuf, "Ambiguous projection specified: %s or %s?", theproj->name, i->name);
+ *error = errbuf;
+ return;
+ }
+ if(*n != i->npar) {
+ sprintf(errbuf, "%s projection requires %d parameter%s", i->name, i->npar, i->npar>1?"s":"");
+ *error = errbuf;
+ return;
+ }
+ if(strcmp(i->name, "bicentric") == 0 ||
+ strcmp(i->name, "elliptic") == 0)
+ par[0] = -par[0];
+ switch(*n) {
+ case 0: projfun = (i->prog)(); break;
+ case 1: projfun = (i->prog)(par[0]); break;
+ case 2: projfun = (i->prog)(par[0], par[1]); break;
+ }
+ theproj = i;
+ }
+ }
+ if(theproj == 0) {
+ sprintf(errbuf, "Unknown projection: %s", *name);
+ *error = errbuf;
+ return;
+ }
+ orient(o[0], -o[1], -o[2]);
+}
+
+static int project(lon, lat, x, y)
+ double lon, lat;
+ double *x, *y;
+{
+ struct place p;
+
+ if(projfun == 0) {
+ *x = lon;
+ *y = lat;
+ return(1);
+ }
+ lon = DEG2RAD(lon);
+ lat = DEG2RAD(lat);
+ p.wlon.l = -lon; p.wlon.s = sin(-lon); p.wlon.c = cos(-lon);
+ p.nlat.l = lat; p.nlat.s = sin(lat); p.nlat.c = cos(lat);
+ normalize(&p);
+ return((*projfun)(&p, x, y));
+}
+
+void doproj(lon, lat, n, range, error)
+ double lon[], lat[], range[];
+ long *n, *error;
+{
+ int i, ok;
+ double x, y;
+
+ *error = 0;
+ range[XMIN] = range[YMIN] = SINGLE_XMAX;
+ range[XMAX] = range[YMAX] = -SINGLE_XMAX;
+ for(i = 0; i < *n; i++, lon++, lat++) {
+ if(ISNA(*lon) || ISNA(*lat))
+ continue;
+ ok = 1 == project(*lon, *lat, &x, &y);
+ if(!ok || ABS(x) > SINGLE_XMAX || ABS(y) > SINGLE_XMAX) {
+ *error = 1;
+ *lon = NA_REAL;
+ *lat = NA_REAL;
+ } else {
+ *lon = x;
+ *lat = y;
+ range[XMIN] = MIN(range[XMIN], x);
+ range[XMAX] = MAX(range[XMAX], x);
+ range[YMIN] = MIN(range[YMIN], y);
+ range[YMAX] = MAX(range[YMAX], y);
+ }
+ }
+}
+
diff --git a/src/mercator.c b/src/mercator.c
new file mode 100644
index 0000000..13aa787
--- /dev/null
+++ b/src/mercator.c
@@ -0,0 +1,41 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static int
+Xmercator(struct place *place, double *x, double *y)
+{
+ if(fabs(place->nlat.l) > 80.*RAD)
+ return(-1);
+ *x = -place->wlon.l;
+ *y = 0.5*log((1+place->nlat.s)/(1-place->nlat.s));
+ return(1);
+}
+
+proj
+mercator(void)
+{
+ return(Xmercator);
+}
+
+static double ecc = ECC;
+
+static int
+Xspmercator(struct place *place, double *x, double *y)
+{
+ if(Xmercator(place,x,y) < 0)
+ return(-1);
+ *y += 0.5*ecc*log((1-ecc*place->nlat.s)/(1+ecc*place->nlat.s));
+ return(1);
+}
+
+proj
+sp_mercator(void)
+{
+ return(Xspmercator);
+}
diff --git a/src/mollweide.c b/src/mollweide.c
new file mode 100644
index 0000000..05b9fca
--- /dev/null
+++ b/src/mollweide.c
@@ -0,0 +1,30 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static int
+Xmollweide(struct place *place, double *x, double *y)
+{
+ double z;
+ double w;
+ z = place->nlat.l;
+ if(fabs(z)<89.9*RAD)
+ do { /*newton for 2z+sin2z=pi*sin(lat)*/
+ w = (2*z+sin(2*z)-PI*place->nlat.s)/(2+2*cos(2*z));
+ z -= w;
+ } while(fabs(w)>=.00001);
+ *y = sin(z);
+ *x = - (2/PI)*cos(z)*place->wlon.l;
+ return(1);
+}
+
+proj
+mollweide(void)
+{
+ return(Xmollweide);
+}
diff --git a/src/newyorker.c b/src/newyorker.c
new file mode 100644
index 0000000..06072a7
--- /dev/null
+++ b/src/newyorker.c
@@ -0,0 +1,33 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static double a;
+
+static int
+Xnewyorker(struct place *place, double *x, double *y)
+{
+ double r = PI/2 - place->nlat.l;
+ double s;
+ if(r<.001) /* cheat to plot center */
+ s = 0;
+ else if(r<a)
+ return -1;
+ else
+ s = log(r/a);
+ *x = -s * place->wlon.s;
+ *y = -s * place->wlon.c;
+ return(1);
+}
+
+proj
+newyorker(double a0)
+{
+ a = a0*RAD;
+ return(Xnewyorker);
+}
diff --git a/src/orthographic.c b/src/orthographic.c
new file mode 100644
index 0000000..64e38dd
--- /dev/null
+++ b/src/orthographic.c
@@ -0,0 +1,40 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+
+int
+Xorthographic(struct place *place, double *x, double *y)
+{
+ *x = - place->nlat.c * place->wlon.s;
+ *y = - place->nlat.c * place->wlon.c;
+ return(place->nlat.l<0.? 0 : 1);
+}
+
+proj
+orthographic(void)
+{
+ return(Xorthographic);
+}
+
+int
+olimb(double *lat, double *lon, double res)
+{
+ static int first = 1;
+ if(first) {
+ *lat = 0;
+ *lon = -180;
+ first = 0;
+ return 0;
+ }
+ *lon += res;
+ if(*lon <= 180)
+ return 1;
+ first = 1;
+ return -1;
+}
diff --git a/src/perspective.c b/src/perspective.c
new file mode 100644
index 0000000..382067f
--- /dev/null
+++ b/src/perspective.c
@@ -0,0 +1,89 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+#define ORTHRAD 1000
+static double viewpt;
+
+static int
+Xperspective(struct place *place, double *x, double *y)
+{
+ double r;
+ if(viewpt<=1+FUZZ && fabs(place->nlat.s)<=viewpt+.01)
+ return(-1);
+ r = place->nlat.c*(viewpt - 1.)/(viewpt - place->nlat.s);
+ *x = - r*place->wlon.s;
+ *y = - r*place->wlon.c;
+ if(r>4.)
+ return(-1);
+ if((fabs(viewpt)>1 && place->nlat.s<1/viewpt) ||
+ (fabs(viewpt)<=1 && place->nlat.s<viewpt))
+ return 0;
+ return(1);
+}
+
+proj
+map_perspective(double radius)
+{
+ viewpt = radius;
+ if(viewpt >= ORTHRAD)
+ return(Xorthographic);
+ if(fabs(viewpt-1.)<.0001)
+ return(0);
+ return(Xperspective);
+}
+
+ /* called from various conformal projections,
+ but not from stereographic itself */
+int
+Xstereographic(struct place *place, double *x, double *y)
+{
+ double v = viewpt;
+ int retval;
+ viewpt = -1;
+ retval = Xperspective(place, x, y);
+ viewpt = v;
+ return retval;
+}
+
+proj
+stereographic(void)
+{
+ viewpt = -1.;
+ return(Xperspective);
+}
+
+proj
+gnomonic(void)
+{
+ viewpt = 0.;
+ return(Xperspective);
+}
+
+int
+plimb(double *lat, double *lon, double res)
+{
+ static int first = 1;
+ if(viewpt >= ORTHRAD)
+ return olimb(lat, lon, res);
+ if(first) {
+ first = 0;
+ *lon = -180;
+ if(fabs(viewpt) < .01)
+ *lat = 0;
+ else if(fabs(viewpt)<=1)
+ *lat = asin(viewpt)/RAD;
+ else
+ *lat = asin(1/viewpt)/RAD;
+ } else
+ *lon += res;
+ if(*lon <= 180)
+ return 1;
+ first = 1;
+ return -1;
+}
diff --git a/src/polyconic.c b/src/polyconic.c
new file mode 100644
index 0000000..626b5ea
--- /dev/null
+++ b/src/polyconic.c
@@ -0,0 +1,33 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+int
+Xpolyconic(struct place *place, double *x, double *y)
+{
+ double r, alpha;
+ double lat2, lon2;
+ if(fabs(place->nlat.l) > .01) {
+ r = place->nlat.c / place->nlat.s;
+ alpha = place->wlon.l * place->nlat.s;
+ *y = place->nlat.l + r*(1 - cos(alpha));
+ *x = - r*sin(alpha);
+ } else {
+ lon2 = place->wlon.l * place->wlon.l;
+ lat2 = place->nlat.l * place->nlat.l;
+ *y = place->nlat.l * (1+(lon2/2)*(1-(8+lon2)*lat2/12));
+ *x = - place->wlon.l * (1-lat2*(3+lon2)/6);
+ }
+ return(1);
+}
+
+proj
+polyconic(void)
+{
+ return(Xpolyconic);
+}
diff --git a/src/rectangular.c b/src/rectangular.c
new file mode 100644
index 0000000..8367956
--- /dev/null
+++ b/src/rectangular.c
@@ -0,0 +1,27 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static double scale;
+
+static int
+Xrectangular(struct place *place, double *x, double *y)
+{
+ *x = -scale*place->wlon.l;
+ *y = place->nlat.l;
+ return(1);
+}
+
+proj
+rectangular(double par)
+{
+ scale = cos(par*RAD);
+ if(scale<.1)
+ return 0;
+ return(Xrectangular);
+}
diff --git a/src/simpleconic.c b/src/simpleconic.c
new file mode 100644
index 0000000..46aecb2
--- /dev/null
+++ b/src/simpleconic.c
@@ -0,0 +1,39 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static double r0, a;
+
+static int
+Xsimpleconic(struct place *place, double *x, double *y)
+{
+ double r = r0 - place->nlat.l;
+ double t = a*place->wlon.l;
+ *x = -r*sin(t);
+ *y = -r*cos(t);
+ return 1;
+}
+
+proj
+simpleconic(double par0, double par1)
+{
+ struct coord lat0;
+ struct coord lat1;
+ deg2rad(par0,&lat0);
+ deg2rad(par1,&lat1);
+ if(fabs(lat0.l+lat1.l)<.01)
+ return rectangular(par0);
+ if(fabs(lat0.l-lat1.l)<.01) {
+ a = lat0.s/lat0.l;
+ r0 = lat0.c/lat0.s + lat0.l;
+ } else {
+ a = (lat1.c-lat0.c)/(lat0.l-lat1.l);
+ r0 = ((lat0.c+lat1.c)/a + lat1.l + lat0.l)/2;
+ }
+ return Xsimpleconic;
+}
diff --git a/src/sinusoidal.c b/src/sinusoidal.c
new file mode 100644
index 0000000..5b67c3e
--- /dev/null
+++ b/src/sinusoidal.c
@@ -0,0 +1,22 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+int
+Xsinusoidal(struct place *place, double *x, double *y)
+{
+ *x = - place->wlon.l * place->nlat.c;
+ *y = place->nlat.l;
+ return(1);
+}
+
+proj
+sinusoidal(void)
+{
+ return(Xsinusoidal);
+}
diff --git a/src/tetra.c b/src/tetra.c
new file mode 100644
index 0000000..c15230b
--- /dev/null
+++ b/src/tetra.c
@@ -0,0 +1,212 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+/*
+ * conformal map of earth onto tetrahedron
+ * the stages of mapping are
+ * (a) stereo projection of tetrahedral face onto
+ * isosceles curvilinear triangle with 3 120-degree
+ * angles and one straight side
+ * (b) map of this triangle onto half plane cut along
+ * 3 rays from the roots of unity to infinity
+ * formula (z^4+2*3^.5*z^2-1)/(z^4-2*3^.5*z^2-1)
+ * (c) do 3 times for each sector of plane:
+ * map of |arg z|<=pi/6, cut along z>1 into
+ * triangle |arg z|<=pi/6, Re z<=const,
+ * with upper side of cut going into upper half of
+ * of vertical side of triangle and lowere into lower
+ * formula int from 0 to z dz/sqrt(1-z^3)
+ *
+ * int from u to 1 3^.25*du/sqrt(1-u^3) =
+ F(acos((rt3-1+u)/(rt3+1-u)),sqrt(1/2+rt3/4))
+ * int from 1 to u 3^.25*du/sqrt(u^3-1) =
+ * F(acos((rt3+1-u)/(rt3-1+u)),sqrt(1/2-rt3/4))
+ * this latter formula extends analytically down to
+ * u=0 and is the basis of this routine, with the
+ * argument of complex elliptic integral elco2
+ * being tan(acos...)
+ * the formula F(pi-x,k) = 2*F(pi/2,k)-F(x,k) is
+ * used to cross over into the region where Re(acos...)>pi/2
+ * f0 and fpi are suitably scaled complete integrals
+*/
+
+#define TFUZZ 0.00001
+
+static struct place tpole[4]; /* point of tangency of tetrahedron face*/
+static double tpoleinit[4][2] = {
+ { 1., 0. },
+ { 1., 180.},
+ { -1., 90.},
+ { -1., -90.}
+};
+static struct tproj {
+ double tlat,tlon; /* center of stereo projection*/
+ double ttwist; /* rotatn before stereo*/
+ double trot; /*rotate after projection*/
+ struct place projpl; /*same as tlat,tlon*/
+ struct coord projtw; /*same as ttwist*/
+ struct coord postrot; /*same as trot*/
+} tproj[4][4] = {
+{/*00*/ {0.},
+ /*01*/ {90., 0., 90., -90.},
+ /*02*/ {0., 45., -45., 150.},
+ /*03*/ {0., -45., -135., 30.}
+},
+{/*10*/ {90., 0., -90., 90.},
+ /*11*/ {0.},
+ /*12*/ {0., 135., -135., -150.},
+ /*13*/ {0., -135., -45., -30.}
+},
+{/*20*/ {0., 45., 135., -30.},
+ /*21*/ {0., 135., 45., -150.},
+ /*22*/ {0.},
+ /*23*/ {-90., 0., 180., 90.}
+},
+{/*30*/ {0., -45., 45., -150.},
+ /*31*/ {0., -135., 135., -30.},
+ /*32*/ {-90., 0., 0., 90.},
+ /*33*/ {0.}
+}};
+static double tx0[4] = { /*where to move facet after final rotation*/
+ 0., 0., -1., 1. /*-1,1 to be sqrt(3)*/
+};
+static double ty0[4] = {
+ 0., 2., -1., -1.
+};
+static double tx[4];
+static double ty[4];
+static double root3;
+static double rt3inv;
+static double two_rt3;
+static double tkc,tk,tcon;
+static double f0r,f0i,fpir,fpii;
+
+static void
+twhichp(struct place *g, int *p, int *q)
+{
+ int i,j,k;
+ double cosdist[4];
+ struct place *tp;
+ for(i=0;i<4;i++) {
+ tp = &tpole[i];
+ cosdist[i] = g->nlat.s*tp->nlat.s +
+ g->nlat.c*tp->nlat.c*(
+ g->wlon.s*tp->wlon.s +
+ g->wlon.c*tp->wlon.c);
+ }
+ j = 0;
+ for(i=1;i<4;i++)
+ if(cosdist[i] > cosdist[j])
+ j = i;
+ *p = j;
+ k = j==0?1:0;
+ for(i=0;i<4;i++)
+ if(i!=j&&cosdist[i]>cosdist[k])
+ k = i;
+ *q = k;
+}
+
+int
+Xtetra(struct place *place, double *x, double *y)
+{
+ int i,j;
+ struct place pl;
+ register struct tproj *tpp;
+ double vr, vi;
+ double br, bi;
+ double zr,zi,z2r,z2i,z4r,z4i,sr,si,tr,ti;
+ twhichp(place,&i,&j);
+ copyplace(place,&pl);
+ norm(&pl,&tproj[i][j].projpl,&tproj[i][j].projtw);
+ Xstereographic(&pl,&vr,&vi);
+ zr = vr/2;
+ zi = vi/2;
+ if(zr<=TFUZZ)
+ zr = TFUZZ;
+ csq(zr,zi,&z2r,&z2i);
+ csq(z2r,z2i,&z4r,&z4i);
+ z2r *= two_rt3;
+ z2i *= two_rt3;
+ cdiv(z4r+z2r-1,z4i+z2i,z4r-z2r-1,z4i-z2i,&sr,&si);
+ map_csqrt(sr-1,si,&tr,&ti);
+ cdiv(tcon*tr,tcon*ti,root3+1-sr,-si,&br,&bi);
+ if(br<0) {
+ br = -br;
+ bi = -bi;
+ if(!elco2(br,bi,tk,1.,1.,&vr,&vi))
+ return 0;
+ vr = fpir - vr;
+ vi = fpii - vi;
+ } else
+ if(!elco2(br,bi,tk,1.,1.,&vr,&vi))
+ return 0;
+ if(si>=0) {
+ tr = f0r - vi;
+ ti = f0i + vr;
+ } else {
+ tr = f0r + vi;
+ ti = f0i - vr;
+ }
+ tpp = &tproj[i][j];
+ *x = tr*tpp->postrot.c +
+ ti*tpp->postrot.s + tx[i];
+ *y = ti*tpp->postrot.c -
+ tr*tpp->postrot.s + ty[i];
+ return(1);
+}
+
+int
+tetracut(struct place *g, struct place *og, double *cutlon)
+{
+ int i,j,k;
+ if((g->nlat.s<=-rt3inv&&og->nlat.s<=-rt3inv) &&
+ (ckcut(g,og,*cutlon=0.)==2||ckcut(g,og,*cutlon=PI)==2))
+ return(2);
+ twhichp(g,&i,&k);
+ twhichp(og,&j,&k);
+ if(i==j||i==0||j==0)
+ return(1);
+ return(0);
+}
+
+proj
+tetra(void)
+{
+ int i,j;
+ register struct place *tp;
+ register struct tproj *tpp;
+ double t;
+ root3 = sqrt(3.);
+ rt3inv = 1/root3;
+ two_rt3 = 2*root3;
+ tkc = sqrt(.5-.25*root3);
+ tk = sqrt(.5+.25*root3);
+ tcon = 2*sqrt(root3);
+ elco2(tcon/(root3-1),0.,tkc,1.,1.,&f0r,&f0i);
+ elco2(1.e15,0.,tk,1.,1.,&fpir,&fpii);
+ fpir *= 2;
+ fpii *= 2;
+ for(i=0;i<4;i++) {
+ tx[i] = tx0[i]*f0r*root3;
+ ty[i] = ty0[i]*f0r;
+ tp = &tpole[i];
+ t = tp->nlat.s = tpoleinit[i][0]/root3;
+ tp->nlat.c = sqrt(1 - t*t);
+ tp->nlat.l = atan2(tp->nlat.s,tp->nlat.c);
+ deg2rad(tpoleinit[i][1],&tp->wlon);
+ for(j=0;j<4;j++) {
+ tpp = &tproj[i][j];
+ latlon(tpp->tlat,tpp->tlon,&tpp->projpl);
+ deg2rad(tpp->ttwist,&tpp->projtw);
+ deg2rad(tpp->trot,&tpp->postrot);
+ }
+ }
+ return(Xtetra);
+}
+
diff --git a/src/trapezoidal.c b/src/trapezoidal.c
new file mode 100644
index 0000000..e3d4693
--- /dev/null
+++ b/src/trapezoidal.c
@@ -0,0 +1,35 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static struct coord stdpar0, stdpar1;
+static double k;
+static double yeq;
+
+static int
+Xtrapezoidal(struct place *place, double *x, double *y)
+{
+ *y = yeq + place->nlat.l;
+ *x = *y*k*place->wlon.l;
+ return 1;
+}
+
+proj
+trapezoidal(double par0, double par1)
+{
+ if(fabs(fabs(par0)-fabs(par1))<.1)
+ return rectangular(par0);
+ deg2rad(par0,&stdpar0);
+ deg2rad(par1,&stdpar1);
+ if(fabs(par1-par0) < .1)
+ k = stdpar1.s;
+ else
+ k = (stdpar1.c-stdpar0.c)/(stdpar0.l-stdpar1.l);
+ yeq = -stdpar1.l - stdpar1.c/k;
+ return Xtrapezoidal;
+}
diff --git a/src/twocirc.c b/src/twocirc.c
new file mode 100644
index 0000000..733cc17
--- /dev/null
+++ b/src/twocirc.c
@@ -0,0 +1,85 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static double
+quadratic(double a, double b, double c)
+{
+ double disc = b*b - 4*a*c;
+ return disc<0? 0: (-b-sqrt(disc))/(2*a);
+}
+
+/* for projections with meridians being circles centered
+on the x axis and parallels being circles centered on the y
+axis. Find the intersection of the meridian thru (m,0), (90,0),
+with the parallel thru (0,p), (p1,p2) */
+
+static int
+twocircles(double m, double p, double p1, double p2, double *x, double *y)
+{
+ double a; /* center of meridian circle, a>0 */
+ double b; /* center of parallel circle, b>0 */
+ double t,bb;
+ if(m > 0) {
+ twocircles(-m,p,p1,p2,x,y);
+ *x = -*x;
+ } else if(p < 0) {
+ twocircles(m,-p,p1,-p2,x,y);
+ *y = -*y;
+ } else if(p < .01) {
+ *x = m;
+ t = m/p1;
+ *y = p + (p2-p)*t*t;
+ } else if(m > -.01) {
+ *y = p;
+ *x = m - m*p*p;
+ } else {
+ b = p>=1? 1: p>.99? 0.5*(p+1 + p1*p1/(1-p)):
+ 0.5*(p*p-p1*p1-p2*p2)/(p-p2);
+ a = .5*(m - 1/m);
+ t = m*m-p*p+2*(b*p-a*m);
+ bb = b*b;
+ *x = quadratic(1+a*a/bb, -2*a + a*t/bb,
+ t*t/(4*bb) - m*m + 2*a*m);
+ *y = (*x*a+t/2)/b;
+ }
+ return 1;
+}
+
+static int
+Xglobular(struct place *place, double *x, double *y)
+{
+ twocircles(-2*place->wlon.l/PI,
+ 2*place->nlat.l/PI, place->nlat.c, place->nlat.s, x, y);
+ return 1;
+}
+
+proj
+globular(void)
+{
+ return Xglobular;
+}
+
+static int
+Xvandergrinten(struct place *place, double *x, double *y)
+{
+ double t = 2*place->nlat.l/PI;
+ double abst = fabs(t);
+ double pval = abst>=1? 1: abst/(1+sqrt(1-t*t));
+ double p2 = 2*pval/(1+pval);
+ twocircles(-place->wlon.l/PI, pval, sqrt(1-p2*p2), p2, x, y);
+ if(t < 0)
+ *y = -*y;
+ return 1;
+}
+
+proj
+vandergrinten(void)
+{
+ return Xvandergrinten;
+}
diff --git a/src/zcoord.c b/src/zcoord.c
new file mode 100644
index 0000000..a124f0c
--- /dev/null
+++ b/src/zcoord.c
@@ -0,0 +1,149 @@
+/************************************************************
+
+Copyright (C) 1998, Lucent Technologies
+All rights reserved
+
+************************************************************/
+
+#include "map.h"
+
+static double cirmod(double);
+
+static struct place pole; /* map pole is tilted to here */
+static struct coord twist; /* then twisted this much */
+static struct place ipole; /* inverse transfrom */
+static struct coord itwist;
+
+void
+orient(double lat, double lon, double theta)
+{
+ lat = cirmod(lat);
+ if(lat>90.) {
+ lat = 180. - lat;
+ lon -= 180.;
+ theta -= 180.;
+ } else if(lat < -90.) {
+ lat = -180. - lat;
+ lon -= 180.;
+ theta -= 180;
+ }
+ latlon(lat,lon,&pole);
+ deg2rad(theta, &twist);
+ latlon(lat,180.-theta,&ipole);
+ deg2rad(180.-lon, &itwist);
+}
+
+void
+latlon(double lat, double lon, struct place *p)
+{
+ lat = cirmod(lat);
+ if(lat>90.) {
+ lat = 180. - lat;
+ lon -= 180.;
+ } else if(lat < -90.) {
+ lat = -180. - lat;
+ lon -= 180.;
+ }
+ deg2rad(lat,&p->nlat);
+ deg2rad(lon,&p->wlon);
+}
+
+void
+deg2rad(double theta, struct coord *coord)
+{
+ theta = cirmod(theta);
+ coord->l = theta*RAD;
+ if(theta==90) {
+ coord->s = 1;
+ coord->c = 0;
+ } else if(theta== -90) {
+ coord->s = -1;
+ coord->c = 0;
+ } else
+ trig(coord);
+}
+
+static double
+cirmod(double theta)
+{
+ while(theta >= 180.)
+ theta -= 360;
+ while(theta<-180.)
+ theta += 360.;
+ return(theta);
+}
+
+double
+trigclamp(double x)
+{
+ return x>1? 1: x<-1? -1: x;
+}
+
+void
+trig(struct coord *coord)
+{
+ coord->s = sin(coord->l);
+ coord->c = cos(coord->l);
+}
+
+void
+normalize(struct place *gg)
+{
+ norm(gg,&pole,&twist);
+}
+
+void
+invert(struct place *g)
+{
+ norm(g,&ipole,&itwist);
+}
+
+void
+norm(struct place *gg, struct place *pp, struct coord *tw)
+{
+ register struct place *g; /*geographic coords */
+ register struct place *p; /* new pole in old coords*/
+ struct place m; /* standard map coords*/
+ g = gg;
+ p = pp;
+ if(p->nlat.s == 1.) {
+ if(p->wlon.l+tw->l == 0.)
+ return;
+ g->wlon.l -= p->wlon.l+tw->l;
+ } else {
+ if(p->wlon.l != 0) {
+ g->wlon.l -= p->wlon.l;
+ trig(&g->wlon);
+ }
+ m.nlat.s = trigclamp(
+ p->nlat.s * g->nlat.s
+ + p->nlat.c * g->nlat.c * g->wlon.c);
+ m.nlat.c = sqrt(1. - m.nlat.s * m.nlat.s);
+ m.nlat.l = atan2(m.nlat.s, m.nlat.c);
+ m.wlon.s = g->nlat.c * g->wlon.s;
+ m.wlon.c = trigclamp(
+ p->nlat.c * g->nlat.s
+ - p->nlat.s * g->nlat.c * g->wlon.c);
+ m.wlon.l = atan2(m.wlon.s, - m.wlon.c)
+ - tw->l;
+ *g = m;
+ }
+ trig(&g->wlon);
+ if(g->wlon.l>PI)
+ g->wlon.l -= 2*PI;
+ else if(g->wlon.l<-PI)
+ g->wlon.l += 2*PI;
+}
+
+void
+printp(struct place *g)
+{
+printf("%.3f %.3f %.3f %.3f %.3f %.3f\n",
+g->nlat.l,g->nlat.s,g->nlat.c,g->wlon.l,g->wlon.s,g->wlon.c);
+}
+
+void
+copyplace(struct place *g1, struct place *g2)
+{
+ *g2 = *g1;
+}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-mapproj.git
More information about the debian-science-commits
mailing list