void pj_clear_initcache()
{
  if( cache_alloc > 0 )
  {
    int i;

    pj_acquire_lock();

    for( i = 0; i < cache_count; i++ )
      {
	paralist *n, *t = cache_paralist[i];
		
	pj_dalloc( cache_key[i] );

	/* free parameter list elements */
	for (; t != NULL; t = n) {
	  n = t->next;
	  pj_dalloc(t);
	}
      }

    pj_dalloc( cache_key );
    pj_dalloc( cache_paralist );
    cache_count = 0;
    cache_alloc= 0;
    cache_key = NULL;
    cache_paralist = NULL;

    pj_release_lock();
  }
}

	void /* free 2D array */
freev2(void **v, int nrows) {
	if (v) {
		for (v += nrows; nrows > 0; --nrows)
			pj_dalloc(*--v);
		pj_dalloc(v);
	}
}

void nad_free(struct CTABLE *ct) 
{
    if (ct) {
        if( ct->cvs != NULL )
            pj_dalloc(ct->cvs);

        pj_dalloc(ct);
    }
}

	int
bchgen(projUV a, projUV b, int nu, int nv, projUV **f, projUV(*func)(projUV)) {
	int i, j, k;
	projUV arg, *t, bma, bpa, *c;
	double d, fac;

	bma.u = 0.5 * (b.u - a.u); bma.v = 0.5 * (b.v - a.v);
	bpa.u = 0.5 * (b.u + a.u); bpa.v = 0.5 * (b.v + a.v);
	for ( i = 0; i < nu; ++i) {
		arg.u = cos(PI * (i + 0.5) / nu) * bma.u + bpa.u;
		for ( j = 0; j < nv; ++j) {
			arg.v = cos(PI * (j + 0.5) / nv) * bma.v + bpa.v;
			f[i][j] = (*func)(arg);
			if ((f[i][j]).u == HUGE_VAL)
				return(1);
		}
	}
	if (!(c = (projUV *) vector1(nu, sizeof(projUV)))) return 1;
	fac = 2. / nu;
	for ( j = 0; j < nv ; ++j) {
		for ( i = 0; i < nu; ++i) {
			arg.u = arg.v = 0.;
			for (k = 0; k < nu; ++k) {
				d = cos(PI * i * (k + .5) / nu);
				arg.u += f[k][j].u * d;
				arg.v += f[k][j].v * d;
			}
			arg.u *= fac;
			arg.v *= fac;
			c[i] = arg;
		}
		for (i = 0; i < nu; ++i)
			f[i][j] = c[i];
	}
	pj_dalloc(c);
	if (!(c = (projUV*) vector1(nv, sizeof(projUV)))) return 1;
	fac = 2. / nv;
	for ( i = 0; i < nu; ++i) {
		t = f[i];
		for (j = 0; j < nv; ++j) {
			arg.u = arg.v = 0.;
			for (k = 0; k < nv; ++k) {
				d = cos(PI * j * (k + .5) / nv);
				arg.u += t[k].u * d;
				arg.v += t[k].v * d;
			}
			arg.u *= fac;
			arg.v *= fac;
			c[j] = arg;
		}
		f[i] = c;
		c = t;
	}
	pj_dalloc(c);
	return(0);
}

int nad_ctable_load( struct CTABLE *ct, FILE *fid )

{
    int  a_size;

    fseek( fid, sizeof(struct CTABLE), SEEK_SET );

    /* read all the actual shift values */
    a_size = ct->lim.lam * ct->lim.phi;
    ct->cvs = (FLP *) pj_malloc(sizeof(FLP) * a_size);
    if( ct->cvs == NULL 
        || fread(ct->cvs, sizeof(FLP), a_size, fid) != a_size )
    {
        pj_dalloc( ct->cvs );
        ct->cvs = NULL;

        if( getenv("PROJ_DEBUG") != NULL )
        {
            fprintf( stderr, 
            "ctable loading failed on fread() - binary incompatible?\n" );
        }

        pj_errno = -38;
        return 0;
    }

    return 1;
} 

	static void /* convert columns to power series */
cols(projUV **c, projUV **d, int nu, int nv) {
	projUV *sv, **dd;
	int j, k;

	dd = (projUV **)vector2(nu, nv, sizeof(projUV));
	sv = (projUV *)vector1(nv, sizeof(projUV));
	bclear(d, nu, nv);
	bclear(dd, nu, nv);
	bmove(d[0], c[nu-1], nv);
	for (j = nu-2; j >= 1; --j) {
		for (k = nu-j; k >= 1; --k) {
			bmove(sv, d[k], nv);
			submop(d[k], 2., d[k-1], dd[k], nv);
			bmove(dd[k], sv, nv);
		}
		bmove(sv, d[0], nv);
		subop(d[0], c[j], dd[0], nv);
		bmove(dd[0], sv, nv);
	}
	for (j = nu-1; j >= 1; --j)
		subop(d[j], d[j-1], dd[j], nv);
	submop(d[0], .5, c[0], dd[0], nv);
	freev2((void **) dd, nu);
	pj_dalloc(sv);
}

int nad_ctable2_load( projCtx ctx, struct CTABLE *ct, FILE *fid )

{
    int  a_size;

    fseek( fid, 160, SEEK_SET );

    /* read all the actual shift values */
    a_size = ct->lim.lam * ct->lim.phi;
    ct->cvs = (FLP *) pj_malloc(sizeof(FLP) * a_size);
    if( ct->cvs == NULL 
        || fread(ct->cvs, sizeof(FLP), a_size, fid) != a_size )
    {
        pj_dalloc( ct->cvs );
        ct->cvs = NULL;

        if( getenv("PROJ_DEBUG") != NULL )
        {
            fprintf( stderr,
            "ctable2 loading failed on fread() - binary incompatible?\n" );
        }

        pj_ctx_set_errno( ctx, -38 );
        return 0;
    }

    if( !IS_LSB )
    {
        swap_words( ct->cvs, 4, a_size * 2 );
    }

    return 1;
} 

void pj_gridinfo_free( projCtx ctx, PJ_GRIDINFO *gi )

{
    if( gi == NULL )
        return;

    if( gi->child != NULL )
    {
        PJ_GRIDINFO *child, *next;

        for( child = gi->child; child != NULL; child=next)
        {
            next=child->next;
            pj_gridinfo_free( ctx, child );
        }
    }

    if( gi->ct != NULL )
        nad_free( gi->ct );

    free( gi->gridname );
    if( gi->filename != NULL )
        free( gi->filename );

    pj_dalloc( gi );
}

int pj_apply_gridshift( projCtx ctx, const char *nadgrids, int inverse, 
                        long point_count, int point_offset,
                        double *x, double *y, double *z )

{
    PJ_GRIDINFO **gridlist;
    int           grid_count;
    int           ret;
    
    gridlist = pj_gridlist_from_nadgrids( ctx, nadgrids, &grid_count );

    if( gridlist == NULL || grid_count == 0 )
        return ctx->last_errno;

    ret = pj_apply_gridshift_3( ctx, gridlist, grid_count, inverse, 
                                point_count, point_offset, x, y, z );

    /* 
    ** Note this frees the array of grid list pointers, but not the grids
    ** which is as intended.  The grids themselves live on.
    */
    pj_dalloc( gridlist );

    return ret;
}

	static void /* convert row to pover series */
rows(projUV *c, projUV *d, int n) {
	projUV sv, *dd;
	int j, k;

	dd = (projUV *)vector1(n-1, sizeof(projUV));
	sv.u = sv.v = 0.;
	for (j = 0; j < n; ++j) d[j] = dd[j] = sv;
	d[0] = c[n-1];
	for (j = n-2; j >= 1; --j) {
		for (k = n-j; k >= 1; --k) {
			sv = d[k];
			d[k].u = 2. * d[k-1].u - dd[k].u;
			d[k].v = 2. * d[k-1].v - dd[k].v;
			dd[k] = sv;
		}
		sv = d[0];
		d[0].u = -dd[0].u + c[j].u;
		d[0].v = -dd[0].v + c[j].v;
		dd[0] = sv;
	}
	for (j = n-1; j >= 1; --j) {
		d[j].u = d[j-1].u - dd[j].u;
		d[j].v = d[j-1].v - dd[j].v;
	}
	d[0].u = -dd[0].u + .5 * c[0].u;
	d[0].v = -dd[0].v + .5 * c[0].v;
	pj_dalloc(dd);
}

static VALUE method_proj4_initialize_copy(VALUE self, VALUE orig)
{
    RGeo_Proj4Data* self_data;
    projPJ pj;
    RGeo_Proj4Data* orig_data;
    char* str;

    // Clear out any existing value
    self_data = RGEO_PROJ4_DATA_PTR(self);
    pj = self_data->pj;
    if (pj) {
        pj_free(pj);
        self_data->pj = NULL;
        self_data->original_str = Qnil;
    }

    // Copy value from orig
    orig_data = RGEO_PROJ4_DATA_PTR(orig);
    if (!NIL_P(orig_data->original_str)) {
        self_data->pj = pj_init_plus(RSTRING_PTR(orig_data->original_str));
    }
    else {
        str = pj_get_def(orig_data->pj, 0);
        self_data->pj = pj_init_plus(str);
        pj_dalloc(str);
    }
    self_data->original_str = orig_data->original_str;
    self_data->uses_radians = orig_data->uses_radians;

    return self;
}

	void
geod_set(int argc, char **argv) {
	paralist *start = 0, *curr;
	double es;
	char *name;
	int i;

    /* put arguments into internal linked list */
	if (argc <= 0)
		emess(1, "no arguments in initialization list");
	start = curr = pj_mkparam(argv[0]);
	for (i = 1; i < argc; ++i) {
		curr->next = pj_mkparam(argv[i]);
		curr = curr->next;
	}
	/* set elliptical parameters */
	if (pj_ell_set(pj_get_default_ctx(),start, &geod_a, &es)) emess(1,"ellipse setup failure");
	/* set units */
	if ((name = pj_param(NULL,start, "sunits").s) != NULL) {
		char *s;
                struct PJ_UNITS *unit_list = pj_get_units_ref();
		for (i = 0; (s = unit_list[i].id) && strcmp(name, s) ; ++i) ;
		if (!s)
			emess(1,"%s unknown unit conversion id", name);
		fr_meter = 1. / (to_meter = atof(unit_list[i].to_meter));
	} else
		to_meter = fr_meter = 1.;
	geod_f = es/(1 + sqrt(1 - es));
	geod_ini();
	/* check if line or arc mode */
	if (pj_param(NULL,start, "tlat_1").i) {
		double del_S;
#undef f
		phi1 = pj_param(NULL,start, "rlat_1").f;
		lam1 = pj_param(NULL,start, "rlon_1").f;
		if (pj_param(NULL,start, "tlat_2").i) {
			phi2 = pj_param(NULL,start, "rlat_2").f;
			lam2 = pj_param(NULL,start, "rlon_2").f;
			geod_inv();
			geod_pre();
		} else if ((geod_S = pj_param(NULL,start, "dS").f) != 0.) {
			al12 = pj_param(NULL,start, "rA").f;
			geod_pre();
			geod_for();
		} else emess(1,"incomplete geodesic/arc info");
		if ((n_alpha = pj_param(NULL,start, "in_A").i) > 0) {
			if (!(del_alpha = pj_param(NULL,start, "rdel_A").f))
				emess(1,"del azimuth == 0");
		} else if ((del_S = fabs(pj_param(NULL,start, "ddel_S").f)) != 0.) {
			n_S = (int)(geod_S / del_S + .5);
		} else if ((n_S = pj_param(NULL,start, "in_S").i) <= 0)
			emess(1,"no interval divisor selected");
	}
	/* free up linked list */
	for ( ; start; start = curr) {
		curr = start->next;
		pj_dalloc(start);
	}
}

PJ_INIT_INFO proj_init_info(const char *initname){
/******************************************************************************
    Information about a named init file.

    Maximum length of initname is 64.

    Returns PJ_INIT_INFO struct.

    If the init file is not found all members of the return struct are set
    to the empty string.

    If the init file is found, but the metadata is missing, the value is
    set to "Unknown".
******************************************************************************/
    int file_found;
    char param[80], key[74];
    paralist *start, *next;
    PJ_INIT_INFO ininfo;
    PJ_CONTEXT *ctx = pj_get_default_ctx();

    memset(&ininfo, 0, sizeof(PJ_INIT_INFO));

    file_found = pj_find_file(ctx, initname, ininfo.filename, sizeof(ininfo.filename));
    if (!file_found || strlen(initname) > 64) {
        return ininfo;
    }

    /* The initial memset (0) makes strncpy safe here */
    strncpy (ininfo.name, initname, sizeof(ininfo.name) - 1);
    strcpy(ininfo.origin, "Unknown");
    strcpy(ininfo.version, "Unknown");
    strcpy(ininfo.lastupdate, "Unknown");

    strncpy (key, initname, 64); /* make room for ":metadata\0" at the end */
    key[64] = 0;
    strncat(key, ":metadata", 9);
    strcpy(param, "+init=");
    /* The +strlen(param) avoids a cppcheck false positive warning */
    strncat(param + strlen(param), key, sizeof(param)-1-strlen(param));

    start = pj_mkparam(param);
    pj_expand_init(ctx, start);

    if (pj_param(ctx, start, "tversion").i)
        strncpy(ininfo.version, pj_param(ctx, start, "sversion").s, sizeof(ininfo.version) - 1);

    if (pj_param(ctx, start, "torigin").i)
        strncpy(ininfo.origin, pj_param(ctx, start, "sorigin").s, sizeof(ininfo.origin) - 1);

    if (pj_param(ctx, start, "tlastupdate").i)
        strncpy(ininfo.lastupdate, pj_param(ctx, start, "slastupdate").s, sizeof(ininfo.lastupdate) - 1);

    for ( ; start; start = next) {
        next = start->next;
        pj_dalloc(start);
    }

   return ininfo;
}

    /** Returns the PROJ.4 initialization string suitable for use with
    pj_init_plus() that would produce this coordinate system, but with the
    definition expanded as much as possible (for instance +init= and
    +datum= definitions).
    @param options Unused at this point
    */
    std::string get_def(int options=0) const
    {
        char *pj_def = 0;
        pj_def = pj_get_def(pj, options);

        std::string ret = std::string(pj_def);
        pj_dalloc(pj_def);
        return ret;
    }

void
pj_free(PJ *P) {
    if (P) {
        paralist *t = P->params, *n;

        /* free parameter list elements */
        for (t = P->params; t; t = n) {
            n = t->next;
            pj_dalloc(t);
        }

        /* free array of grid pointers if we have one */
        if( P->gridlist != NULL )
            pj_dalloc( P->gridlist );
        
        /* free projection parameters */
        P->pfree(P);
    }
}

void pj_insert_initcache( const char *filekey, const paralist *list )

{
  pj_acquire_lock();

  /* 
  ** Grow list if required.
  */
  if( cache_count == cache_alloc )
    {
      char **cache_key_new;
      paralist **cache_paralist_new;

      cache_alloc = cache_alloc * 2 + 15;

      cache_key_new = (char **) pj_malloc(sizeof(char*) * cache_alloc);
      memcpy( cache_key, cache_key_new, sizeof(char*) * cache_count);
      pj_dalloc( cache_key );
      cache_key = cache_key_new;

      cache_paralist_new = (paralist **) 
	pj_malloc(sizeof(paralist*) * cache_alloc);
      memcpy( cache_paralist_new, cache_paralist, 
	      sizeof(paralist*) * cache_count );
      pj_dalloc( cache_paralist );
      cache_paralist = cache_paralist_new;
    }

  /*
  ** Duplicate the filekey and paralist, and insert in cache.
  */
  cache_key[cache_count] = (char *) pj_malloc(strlen(filekey)+1);
  strcpy( cache_key[cache_count], filekey );

  cache_paralist[cache_count] = pj_clone_paralist( list );

  cache_count++;

  pj_release_lock();
}

Handle<Value> Projection::GetDefinition(Local<String> name, const AccessorInfo& info)
{
    HandleScope     scope;
    Projection      *self = ObjectWrap::Unwrap<Projection>(info.Holder());
    char            *def;
    Handle<Value>   def_obj;
    
    def = pj_get_def(self->pj, 0);
    def_obj = String::New(def);
    
    pj_dalloc(def);
    return scope.Close(def_obj);
}

static VALUE method_proj4_canonical_str(VALUE self)
{
  VALUE result = Qnil;
  projPJ pj = RGEO_PROJ4_DATA_PTR(self)->pj;
  if (pj) {
    char* str = pj_get_def(pj, 0);
    if (str) {
      result = rb_str_new2(str);
      pj_dalloc(str);
    }
  }
  return result;
}

PJ *
pj_init_plus( const char *definition )

{
#define MAX_ARG 200
    char	*argv[MAX_ARG];
    char	*defn_copy;
    int		argc = 0, i;
    PJ	        *result;
    
    /* make a copy that we can manipulate */
    defn_copy = (char *) pj_malloc( strlen(definition)+1 );
    strcpy( defn_copy, definition );

    /* split into arguments based on '+' and trim white space */

    for( i = 0; defn_copy[i] != '\0'; i++ )
    {
        switch( defn_copy[i] )
        {
          case '+':
            if( i == 0 || defn_copy[i-1] == '\0' )
            {
                if( argc+1 == MAX_ARG )
                {
                    pj_errno = -44;
                    return NULL;
                }
                
                argv[argc++] = defn_copy + i + 1;
            }
            break;

          case ' ':
          case '\t':
          case '\n':
            defn_copy[i] = '\0';
            break;

          default:
            /* do nothing */;
        }
    }

    /* perform actual initialization */
    result = pj_init( argc, argv );

    pj_dalloc( defn_copy );

    return result;
}

static projectionObj* msGetProjectNormalized( const projectionObj* p )
{
  int i;
  char* pszNewProj4Def;
  projectionObj* pnew;

  pnew = (projectionObj*)msSmallMalloc(sizeof(projectionObj));
  msInitProjection(pnew);
  msCopyProjection(pnew, (projectionObj*)p);

  if(p->proj == NULL )
      return pnew;

  /* Normalize definition so that msProjectDiffers() works better */
  pszNewProj4Def = pj_get_def( p->proj, 0 );
  msFreeCharArray(pnew->args, pnew->numargs);
  pnew->args = msStringSplit(pszNewProj4Def,'+', &pnew->numargs);
  for(i = 0; i < pnew->numargs; i++)
  {
      /* Remove trailing space */
      if( strlen(pnew->args[i]) > 0 && pnew->args[i][strlen(pnew->args[i])-1] == ' ' )
          pnew->args[i][strlen(pnew->args[i])-1] = '\0';
      /* Remove spurious no_defs or init= */
      if( strcmp(pnew->args[i], "no_defs") == 0 ||
          strncmp(pnew->args[i], "init=", 5) == 0 )
      {
          if( i < pnew->numargs - 1 )
          {
              msFree(pnew->args[i]);
              memmove(pnew->args + i, pnew->args + i + 1,
                      sizeof(char*) * (pnew->numargs - 1 -i ));
          }
          pnew->numargs --;
          i --;
          continue;
      }
  }
  /* Sort the strings so they can be compared */
  qsort(pnew->args, pnew->numargs, sizeof(char*), msProjectSortString);
  /*{
      fprintf(stderr, "'%s' =\n", pszNewProj4Def);
      for(i = 0; i < p->numargs; i++)
          fprintf(stderr, "'%s' ", p->args[i]);
      fprintf(stderr, "\n");
  }*/
  pj_dalloc(pszNewProj4Def);
  
  return pnew;
}