/* /+ RMNLIB - Library of useful routines for C and FORTRAN programming */
/* * Copyright (C) 1975-2001  Division de Recherche en Prevision Numerique */
/* *                          Environnement Canada */
/* * */
/* * This library is free software; you can redistribute it and/or */
/* * modify it under the terms of the GNU Lesser General Public */
/* * License as published by the Free Software Foundation, */
/* * version 2.1 of the License. */
/* * */
/* * This library is distributed in the hope that it will be useful, */
/* * but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU */
/* * Lesser General Public License for more details. */
/* * */
/* * You should have received a copy of the GNU Lesser General Public */
/* * License along with this library; if not, write to the */
/* * Free Software Foundation, Inc., 59 Temple Place - Suite 330, */
/* * Boston, MA 02111-1307, USA. */
/* +/ */
/* **S/R LLFXY - COMPUTES LATITUDE AND LONGITUDE OF A POINT */
/* **S/R XYFLL - COMPUTES THE GRID CO-ORDINATES OF A POINT */

/* AUTHOR   - J.D. HENDERSON  -  FEB 75 */

/* REVISION 001   C. THIBEAULT  -  NOV 79  DOCUMENTATION */

/* REVISION 002   inVisu Translated to php (from Fortan) */


function r_sign(a,b)
{
  x = (a >= 0 ? a : - a);
  return( b >= 0 ? x : -x);
}


function llfxy(x,y,d60,dgrw,nhem) {

  c_b3 = 90.;
  c_b4 = 180.;
  rdtodg=180./Math.PI;
  re=1.866025*6371000/d60;
  re2=re*re;

  dlat=90.;
  dlon=0.;

  if (x != 0. || y != 0.) {

    if (x == 0.) {
      dlon = r_sign(c_b3, y);
    }
    if (x != 0.) {
      dlon = Math.atan(y / x) * rdtodg;
    }
    if (x < 0.) {
      dlon += r_sign(c_b4, y);
    }


    dlon -= dgrw;
    if (dlon > 180.) {
      dlon += -360.;
    }
    if (dlon < -180.) {
      dlon += 360.;
    }

    r2=x*x + y*y;
    dlat=(re2-r2)/(re2+r2);
    dlat=Math.asin(dlat)*rdtodg;

  }


  if (nhem == 2) {
    dlat = -(dlat);
  }
  if (nhem == 2) {
    dlon = -(dlon);
  }

  return new Array(dlat,dlon);
}

function xyfll(dlat,dlon,d60,dgrw,nhem) {
  dgtord = Math.PI/180.;

  re = 11888445.275 / d60;

  var glon = eval(dlon);
  if (nhem == 2) {
    glon = -dlon;
  }
  var glat = eval(dlat);
  if (nhem == 2) {
    glat = -dlat;
  }

  rlon = dgtord * (glon + dgrw);
  rlat = dgtord * glat;

  sinlat = Math.sin(rlat);
  r = re * Math.sqrt((1. - sinlat) / (sinlat + 1.));
  x = eval(r * Math.cos(rlon));
  y = eval(r * Math.sin(rlon));


  return new Array(x,y);
}


// This returns Lat/Lon of any point in the grid
// FIXME hardcoded
function llfninj(ni,nj) {
  x=-3012.49+ni*5;
  y=-5781.87+nj*5;

  return llfxy(x,y,1000,10,1);
}


function ninjfll(lat,lon) {
  xy=xyfll(lat,lon,1000,10,1);

  x=eval(xy[0]);
  y=eval(xy[1]);
  ni=Math.round((x+3012.49)/5);
  nj=Math.round((y+5781.87)/5);


  return new Array(ni,nj);
}