program polething
	implicit none
c
c	This program calculates the PA dependent term of the zodiacal model
c
	character*46 infile
	character*10 head(5)
	integer*4 i,j
	real*4 frame(360,360),amap(360,360),x,y,eps,epscen,gam,ampl,offset,arg,arg1,xx,yy,yyy
c
	data ampl,offset /8800.,1200./
	data frame /129600*0./
c
	infile ='g:\GegenscheinXStars\c3avg_2003-8rolledoff.grd'
	print *,'reading in average map ',infile
	open(12,file=infile,readonly)
	read(12,1)head
1	format(a10)			
	do j=1,360
		read(12,11)(amap(i,j),i=1,360)
	enddo
	close(12)
c
	do j=1,180
		yy=float(j-180)/2.
		do i=1,360
			y=yy
			xx=float(i-180)/2.
			x=xx
			eps=cosd(x)*cosd(y)
			epscen=acosd(eps)
			if(x.lt.0.)x=360.-x	!correct to rotate's domain
			call rotate(0.,+15.,0.,x,y)
			if(x.gt.180.)x=x-180.	!correct back from rotate's domain
			eps=cosd(x)*cosd(y)
			eps=acosd(eps)
			gam=atan2d(sind(x),tand(y))
c
			arg1=45.*exp(-eps/16.)	!was 24.
			arg=abs(gam)*0.6
			yyy=yy+abs(xx)/6.5
			if(yyy.le.-15.)then
			    arg1=arg1*(0.15+0.85*exp(-6.*(cosd(arg)**4)))+25.*exp(-acosd(cosd(2.*x)*cosd(0.7*(y)))/10.)
			else
				arg1=arg1*(0.15+0.85*exp(-6.*(cosd(arg)**4)))*exp(-(yyy+15.)/10.)
				arg1=arg1+25.*exp(-acosd(cosd(2.*x)*cosd(0.7*(y)))/10.)*exp(-(yyy+15.)/8.)	
			endif
			frame(i,j)=1.5*arg1
			if(epscen.lt.15.)frame(i,j)=0.
		enddo
	enddo
c
	do j=181,360
		yy=float(j-180)/2.
		do i=1,360
			y=yy
			xx=float(i-180)/2.
			x=xx
			eps=cosd(x)*cosd(y)
			epscen=acosd(eps)
			if(x.lt.0.)x=360.-x	!correct to rotate's domain
			call rotate(0.,-21.,0.,x,y)
			if(x.gt.180.)x=x-180.	!correct back from rotate's domain
			eps=cosd(x)*cosd(y)
			eps=acosd(eps)
			gam=atan2d(sind(x),tand(y))
c
			arg1=45.*exp(-eps/16.)	!was 24.
			arg=(180.-abs(gam))*0.6
			yyy=yy-abs(xx)/6.5
			if(yyy.ge.20.)then
			    arg1=arg1*(0.15+0.85*exp(-6.*(cosd(arg)**4)))+25.*exp(-acosd(cosd(2.*x)*cosd(0.7*(y)))/10.)!was *(y-3.)))/10.)
			else
				arg1=arg1*(0.15+0.85*exp(-6.*(cosd(arg)**4)))*exp((yyy-20.)/10.)
				arg1=arg1+25.*exp(-acosd(cosd(2.*x)*cosd(0.7*(y)))/10.)*exp((yyy-20.)/8.)!was *(y-3.)))/10.)	
			endif
			frame(i,j)=arg1
			if(epscen.lt.15.)frame(i,j)=0.
		enddo
	enddo
c
	open(13,file='g:\GegenscheinXStars\polething.grd')
	write(13,13)
13	format('DSAA'/'360 360'/'-90 90'/'-90 90'/'-100 100')
	do j=1,360
c		write(13,11)(amap(i,j)-frame(i,j),i=1,360)
		write(13,11)(frame(i,j),i=1,360)
11		format(20f8.2)				
	enddo
	close(13)
	stop 'done'
	end
c*********************************************************************************	
	subroutine rotate(ALFA,BETA,GAMMA,PHI,RLAT)
	implicit none
c INPUTS:
c	ALFA	real	phase angle of the pole Z(new) in old coordinate system
c	BETA	real	polar angle of the pole Z(new) in old coordinate system
c	GAMMA	real	phase angle of X(new) in coordinate system (2)
c	PHI	real	phase angle in old coordinate system
c	RLAT	real	latitude in old coordinate system (-90<=RLAT<=90)
c OUTPUTS:
c	PHI	real	phase angle in new coordinate system (0<=PHI<360)
c	RLAT	real	latitude in new coordinate system
            real        ALFA
            real        BETA
            real        GAMMA
            real        PHI
            real        RLAT

        double precision        SINB
        double precision        COSB
        double precision        SINL
        double precision        COSL
        double precision        SINP
        double precision        COSP
        double precision        COST

        ! Needed for g77 when using goniod.f
        ! What happens if these are available as intrinsic functions??

        double precision        dsind
        double precision        dcosd
        double precision        dacosd
        double precision        datan2d

        if (abs(RLAT) .gt. 90.0)print *,'Bummer: RLAT = ',RLAT 	!call Say('rotate','E','invalid','latitude')

        PHI  = PHI-ALFA                 ! Rotation (1)

        SINB = dsind(dble(BETA))        ! Rotation (2)
        COSB = dcosd(dble(BETA))
        SINL = dsind(dble(RLAT))
        COSL = dcosd(dble(RLAT))
        SINP = dsind(dble(PHI ))
        COSP = dcosd(dble(PHI ))

        COST = SINL*COSB+COSL*SINB*COSP ! cos(polar angle)

        if (dabs(COST) .ne. 1.0d0) then	! changed abs --> dabs from Paul's, probably doesn't matter...
            SINP = COSL*SINP
            COSP = COSL*COSB*COSP-SINL*SINB

        else                            !Point on new Z-axis; PHI not defined

            !-------
            ! The following statement is inserted for the benefit of the
            ! synoptic map drawing program.
            ! If BETA=0 then SINB=0, COSB=+/-1, i.e. the z-axis stays the same
            ! or is flipped by 180 degrees. The input phase angle for points on
            ! the z-axis will remain the same or is changed to 180-(phase angle)
            ! (just as for points not on the z-axis).

            if (SINB .eq. 0.0d0) COSP = COSB*COSP

        end if

        PHI  = datan2d(SINP,COSP)
        PHI  = PHI-GAMMA                ! Rotation (3)
        PHI  = amod(amod(PHI,360.0)+360.0,360.0)!0<=PHI<360
        RLAT = 90.0d0-dacosd(COST)        ! acosd returns between 0 and 180 deg

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