This is proj.c in view mode; [Download] [Up]
/*
** proj.c
** map projection calculation and generation routines
** for project : MapMaker
** using NeXTStep and mach Unix.
** CPSC 414 Assignment No. 4 Project
**
** Programmed by Bradley Head and Thomas Burkholder
** December 1990
*/
#import "proj.h"
#import "pointdata.h"
#import "appkit/graphics.h"
/********************* THE MAP PROJECTIONS (Forward Sol'n) **********************************/
float correct(float f)
{
if (f>=PI)
return f-(2*PI);
if (f < (-(PI)))
return f+(2*PI);
return f;
}
/*
** Map Projection conversion functions
*/
int Cylindrical(Point *posn, ProjParam *init, Point *pp)
{
float coef;
coef = init->radius*cos(init->phi1);
pp->x = (posn->x - init->lon0)*coef;
pp->y = init->radius*sin(posn->y)/cos(init->phi1);
if (pp->x > (PI*coef)) {
pp->x -= ((2*PI)*coef);
return 1;
}
if (pp->x < ((-(PI))*coef)) {
pp->x += ((2*PI)*coef);
return 1;
}
return 0;
}
int EckertIV(Point *posn, ProjParam *init, Point *pp)
{
float coef;
coef = .42224*init->radius*(1 + cos(posn->y - init->phi1));
pp->x = coef*(posn->x - init->lon0);
pp->y = 1.32650*init->radius*sin(posn->y - init->phi1)/cos(init->phi1);
if (pp->x > (PI*coef)) {
pp->x -= ((2*PI)*coef);
return 1;
}
if (pp->x < ((-(PI))*coef)) {
pp->x += ((2*PI)*coef);
return 1;
}
return 0;
}
int Mercator(Point *posn, ProjParam *init, Point *pp)
{
float coef;
coef = init->radius*0.6; /* 0.6 scales to fit into output window with radius = 1.0 */
pp->x =coef*(posn->x - init->lon0);
pp->y = init->radius*log(tan(.785398 + posn->y/2));
if (pp->x > (PI*coef)) {
pp->x -= ((2*PI)*coef);
return 1;
}
if (pp->x < ((-(PI))*coef)) {
pp->x += ((2*PI)*coef);
return 1;
}
return 0;
}
int Mollweide(Point *posn, ProjParam *init, Point *pp)
{
float coef;
coef = .9003*init->radius*cos(posn->y - init->phi1);
pp->x = (posn->x - init->lon0)*coef;
pp->y = 1.4142*init->radius*sin(posn->y - init->phi1);
if (pp->x > (PI*coef)) {
pp->x -= ((2*PI)*coef);
return 1;
}
if (pp->x < ((-(PI))*coef)) {
pp->x += ((2*PI)*coef);
return 1;
}
return 0;
}
int Ortho(Point *posn, ProjParam *init, Point *pp)
{
float cos_c;
float dl;
dl = posn->x - init->lon0;
pp->x = init->radius*cos(posn->y)*sin(dl);
if ((posn->y != -PI) && (posn->y != PI))
pp->y = init->radius*(cos(init->phi1)*sin(posn->y) -
sin(init->phi1)*cos(posn->y)*cos(dl));
else
pp->y = init->radius*cos(posn->y)*cos(dl);
if (posn->y == PI) /* conditionals for line removal */
pp->y = - pp->y;
cos_c = sin(init->phi1)*sin(posn->y)+cos(init->phi1)*cos(posn->y)*cos(dl);
if (cos_c < 0)
pp->pen = SKIP;
return 0;
}
int Sinusoidal(Point *posn, ProjParam *init, Point *pp)
{
float coef;
float xofs;
coef = init->radius*cos(posn->y - init->phi1);
pp->x = (posn->x - init->lon0)*coef ;
pp->y = init->radius*(posn->y - init->phi1);
if (pp->x > (PI*coef)) { /* clipping conditonals for animation */
pp->x -= ((2*PI)*coef);
return 1;
}
if (pp->x < ((-(PI))*coef)) { /* conditionals for animation */
pp->x += ((2*PI)*coef);
return 1;
}
return 0;
}
int Stereo(Point *posn, ProjParam *init, Point *pp)
{
float dl,k,r;
float Len;
dl = posn->x - init->lon0;
r = init->radius*0.6;
k = 2/(1+sin(init->phi1)*sin(posn->y)+cos(init->phi1)*cos(posn->y)*cos(dl));
pp->x = r*k*cos(posn->y)*sin(dl);
pp->y = r*k*(cos(init->phi1)*sin(posn->y) - sin(init->phi1)*cos(posn->y)*cos(dl));
Len = sqrt((pp->x)*(pp->x) + (pp->y)*(pp->y));
if (Len > init->radius*1.2)
pp->pen = SKIP;
return 0;
}
int None(Point *posn, ProjParam *init, Point *pp)
{
float dl;
float Len;
float coef;
coef = init->radius;
dl = posn->x - init->lon0;
pp->x = coef*dl;
pp->y = coef*(posn->y - init->phi1);
if (pp->x > (PI*coef)) { /* clipping condtionals for animation */
pp->x -= ((2*PI)*coef);
return 1;
}
if (pp->x < ((-(PI))*coef)) {
pp->x += ((2*PI)*coef);
return 1;
}
return 0;
}
/****************************************************/
void DoFrame(PointList *framelist) {
/*
** This function creates two frame lines
** for the projections
** all except stereographic and orthographic
*/
float lat, lon;
Point framept;
framept.pencolour = NX_DKGRAY;
for (lon = -180;lon <= 180;lon += 360) {
for (lat = -75;lat <= 75;lat += 5) {
if (lat == 75)
framept.pen = UP;
else
framept.pen = DOWN;
framept.y = lat *toRADS;
framept.x = lon *toRADS;
addToPointList(framelist ,&framept);
}
}
}
void drawCircleFrame(float radius, PointList *pout) {
/*
** This function generates
** a framing circle
** around stereographic and
** the orthographic
** projections.
*/
Point p;
float step, x;
p.pencolour = NX_DKGRAY;
step = radius/100.0;
for (x = -radius;x <= radius;x += step) { /* to top half of circle */
p.y = sqrt(radius*radius - x*x);
p.x = x;
p.pen = DOWN;
addToPointList(pout,&p);
}
for (x = radius-step;x >= (-radius+step); x -= step) { /* do bottom half */
p.y = -(sqrt(radius*radius - x*x));
p.x = x;
p.pen = DOWN;
addToPointList(pout,&p);
}
p.y = 0.0;
p.x = -radius;
p.pen = UP;
addToPointList(pout,&p);
}
/****************************************************/
void DoGrid(PointList *gridlist) {
/*
** This function generates the graticules
** for the output view
** for all the projections
** it increments at 15 degree intervals
*/
int lat, lon;
Point gridpt ;
gridpt.pencolour = NX_LTGRAY;
gridpt.pen = DOWN;
/*
** Compute the lines of latitude
*/
for (lat = -75;lat <= 75;lat += 15) {
for (lon = -180;lon <= 180; lon += 5)
{
if (lon == 180 )
gridpt.pen = UP;
else
gridpt.pen = DOWN;
gridpt.y = lat *toRADS;
gridpt.x = lon *toRADS;
addToPointList(gridlist,&gridpt);
}
}
/*
** Now compute the Meridians
*/
for (lon = -180;lon < 180; lon += 15) {
for (lat = -75;lat <= 75;lat += 5 )
{
if (lat == 75)
gridpt.pen = UP;
else
gridpt.pen = DOWN;
gridpt.y = lat *toRADS;
gridpt.x = lon *toRADS;
addToPointList(gridlist,&gridpt);
}
}
}
void ComputePoints(int (*proj)(Point *posn, ProjParam *init, Point *pp), ProjParam *init,PointList *inlist,PointList *outlist)
{
/*
** This function computes the output
** points given the input list of (lat,lon) pairs
** it calls the appropriate map projection routine
** determines line clipping in the form of skipping
** unnecessary points (points that won't be drawn to output)
** It places the computed points in an outlist list of points
** that can then be plotted
*/
int i ;
int oldsplit,split;
Point pp, old;
Point *pt;
int a =1;
old.pen = UP;
oldsplit = 1;
for (i = gotoPointInList(inlist,0,&pt);(i);i=gotoNextPointInList(inlist,&pt))
{
pp = *pt;
split = (*proj)(pt,init,&pp);
if ((old.pen == DOWN) && (pp.pen == SKIP))
old.pen = UP;
if (split != oldsplit)
old.pen = UP;
if (a)
a = 0;
else
addToPointList(outlist,&old);
old = pp;
oldsplit = split;
}
if (inlist->quantity) {
pp.pen = UP;
addToPointList(outlist,&pp);
}
}
int convertPoints(PointList *pin, PointList *pout, int gridon, ProjParam *init)
{
/* This function is the externally called function
** It handles whether to compute and draw the
** grid and frames and computes the output points
** by calliing ComputePoints with the appropriate
** map projection.
*/
PointList grid, frame; /* point list for the frame around the projection */
ProjParam frameparam; /* frame drawing initialization paramaters */
frameparam.radius = init->radius; /* get the current radius */
frameparam.lon0 = 0.0; /* fix the central meridian for frame computation */
frameparam.phi1 = init->phi1; /* get the current central line of latitude */
init->lon0 = correct(init->lon0); /* correct the offset for animation */
newPointList(&grid); /* generate a new list for the grid */
newPointList(&frame); /* generate a new list for the frame */
DoFrame(&frame); /* compute the frame */
if (gridon)
DoGrid(&grid); /* if want the grid displayed then compute the grid */
switch (init->proj) /* select the appropriate projection to compute */
{
case CYLINDER:
ComputePoints(&Cylindrical,init,&grid,pout);
ComputePoints(&Cylindrical,&frameparam,&frame,pout);
ComputePoints(&Cylindrical,init,pin,pout);
break;
case ECKERT:
ComputePoints(&EckertIV,init,&grid,pout);
ComputePoints(&EckertIV,&frameparam,&frame,pout);
ComputePoints(&EckertIV,init,pin,pout);
break;
case MERCATOR:
ComputePoints(&Mercator,init,&grid,pout);
ComputePoints(&Mercator,&frameparam,&frame,pout);
ComputePoints(&Mercator,init,pin,pout);
break;
case MOLLWEIDE:
ComputePoints(&Mollweide,init,&grid,pout);
ComputePoints(&Mollweide,&frameparam,&frame,pout);
ComputePoints(&Mollweide,init,pin,pout);
break;
case ORTHO:
ComputePoints(&Ortho,init,&grid,pout);
drawCircleFrame(init->radius,pout);
ComputePoints(&Ortho,init,pin,pout);
break;
case SINUSOIDAL:
ComputePoints(&Sinusoidal,init,&grid,pout);
ComputePoints(&Sinusoidal,&frameparam,&frame,pout);
ComputePoints(&Sinusoidal,init,pin, pout);
break;
case STEREO:
ComputePoints(&Stereo,init,&grid,pout);
drawCircleFrame(init->radius*1.2, pout);
ComputePoints(&Stereo,init,pin,pout);
break;
case NONE:
ComputePoints(&None,init,&grid,pout);
ComputePoints(&None,init,pin,pout);
break;
default: break;
}
return 1;
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.