program dumbbell

	do j=1,180
	    y = j-90.5			! ecl lat
	    do i=-179,180
		x = i-0.5		! sun-centered ecl lng
                f = fnc_dumbbell(x,y)
	        write (*, '(2F6.1,F13.7)') x, y, f
            end do
        end do

        end

        function fnc_dumbbell(x,y)

        ! This program calculates the PA dependent term of the zodiacal model
        ! Rewrite of part of Andy's polething.for

	xx = x
	yy = y

	! Rotate to a coordinate system with the x-axis 15 deg south of
	! the Sun (instead of straight at the Sun)

	beta = 15.0
        if (y .gt. 0.0) beta = -21.0
	call rotate(0.0,beta,0.0,xx,yy)

	eps =  acosd(cosd(xx)*cosd(yy))	! elongation relative to (what?)
	gam = atan2d(sind(xx),tand(yy))	! position angle relative to (what?)
 
        aa  = 1.5                       ! South hemisphere
	gam = abs(gam)
	xs  = abs(x)/6.5-abs(y)+15.0
	if (y .gt. 0.0) then            ! North hemisphere
            aa  = 1.0
	    gam = 180.0-gam
	    xs  = xs+5.0
	end if

	eps = 60.0*exp(-eps/16.0)*(0.15+0.85*exp(-6.0*(cosd(gam*0.6)**4)))
	gam = 25.0*exp(-acosd(cosd(2.0*xx)*cosd(0.7*yy))/10.0)

	if (xs .le. 0.0)then
	    eps = eps+gam
	else
	    eps = eps*exp(-xs/10.0)+gam*exp(-xs/8.0)
	end if

        fnc_dumbbell = aa*eps

        return 
	end