c *************************************************************
       Program lambert
c *****************************************************************
c * This subroutine finds the geographical coordinates of a       *
c * given pixel of a Lambert conformal conic image.               *
c *                                                               *
c * Louis Moreau, Intermap Information Technologies Ltd.          *
c * Updated, April 1996.                                          *
c * Reference: Map Projections - A Working Manual, US Geological  *
c * Survey (paper 1395) by John P. Snyder, 1989.                  *
c *****************************************************************
       implicit none
c       real Nx
        real a,e,phio,lamo,phi,lam,phi1,phi2
        real theta,t
        real m1,m2,t1,t2,to,ro,roo,n,F
        real pi,XMc,YMc,x,y,ksi,tm,Xm,Ym
        real lamc,phic,tc
        character*1 ans
 
100     pi=3.141592654

        write(*,1) 
1       format('Enter the pixel number (1 is the leftmost pixel):')
        read(*,*) x
4       format('Enter the line number (1 is the first line):')
        write(*,4)
        read(*,*) y
        x=x-1
        y=y-1

c ***** GRS 80 ellipsoid ************
c ** values are in km (or pixel assuming 1 pixel/km)****
c a: equatorial axis Earth
c e: excentricity of Earth

       a=6378.137
       e=0.081819221
c *****************************************************
c ** Lambert conformal reference values from GeoComp **
c ** MAY HAVE TO BE CHANGED FOR OTHER LCC PROJECTION **
c phi1, phi2 : standard paralleles
c lamo,phi : reference long. and lat. (deg) i.e. where XM=0 and YM=0

        phi1=49.0
        phi2=77.0
        lamo=-95.
        phio=0.
        phi1=phi1*pi/180.
        phi2=phi2*pi/180.
        lamo=lamo*pi/180.
        phio=phio*pi/180.

c *** Top left corner values i.e. where x=0, y=0******
c ** MAY HAVE TO BE CHANGED FOR OTHER IMAGERIES *****

        lamc=-115.4085917
        phic=59.36395
        phic=phic*pi/180.
        lamc=lamc*pi/180.
c **********************************************************************
c ** XMc,YMc: values of the top left corner 
c ** Note: The XM,YM cartesian system has its origin at phio and lamo
c ** YM increases northward and XM increases eastward. The x,y cartesian
c ** system has its origin at the top left corner of the image (lamc,phic).
c ** y (line) increases southward and x (pixel), eastward.
c ***********************************************************************

        m1=cos(phi1)/((1-e**2*sin(phi1)**2)**0.5)
        m2=cos(phi2)/((1-e**2*sin(phi2)**2)**0.5)
        t1=t(phi1,e)
        t2=t(phi2,e)
        to=t(phio,e)

        n=(log(m1)-log(m2))/(log(t1)-log(t2))
        F=m1/(n*t1**n)           
        roo=a*F*to**n

        tc=t(phic,e)
        ro=a*F*tc**n
        theta=n*(lamc-lamo)
        XMc=ro*sin(theta)
        YMc=roo-ro*cos(theta)
 
c ** x and y are the cartesian coordinates with (0,0) at the
c * top-left corner (in km) **
       
        XM=x - abs(XMc)         
        YM=abs(YMc) - y
c *** find phi, lam **************
c        write(*,*) XM,YM

        ro=(XM**2+(roo-YM)**2)**0.5
        ro=sign(ro,n)
        theta=atan(XM/(roo-YM))
        tm=(ro/a/F)**(1/n)
        ksi=pi/2-2*atan(tm)
        phi=ksi+(e**2/2+5*e**4/24+e**6/12+13*e**8
     + /360)*sin(2*ksi)+(7*e**4/48+29*e**6/240
     + +811*e**8/11520)*sin(4*ksi)+(7*e**6/120+
     + 81*e**8/1120)*sin(6*ksi)+(4279*e**8/161280)
     + *sin(8*ksi)

        lam=theta/n+lamo

        phi=phi*180/pi
        lam=lam*180/pi
        write(*,2) phi
        write(*,3) lam
2       format('Latitude of point (degree):', f9.4)
3       format('Longitude of point (degree):', f9.4)
5       format(//'Another one? (y/n)')
        write(*,5)
        read(*,*) ans
        if((ans.eq.'y').or.(ans.eq.'Y')) goto 100
        end

c **************************************************
        real function t(i,e)
        real i,e

        t=((1-sin(i))/(1+sin(i))*((1+e*sin(i))/
     +      (1-e*sin(i)))**e)**0.5

        return
        end
c ************************************************* 
