PROGRAM genfit
*
      INTEGER NLIM,PMAX
      REAL*8 TINY
      PARAMETER( NLIM=50, PMAX=1024, TINY=1.D-12 )
      COMMON/parm/ ln2,pi,eps,looplim
      REAL*8 ln2,pi,eps
      INTEGER looplim
      COMMON/fac/ fac(0:NLIM)
      REAL*8 fac
      COMMON/prof/ xp(PMAX),yp(PMAX),np
      REAL*8 xp,yp
      INTEGER np
      COMMON/cn/ cn(0:NLIM),cnp(0:NLIM),bnk(0:NLIM,0:NLIM),pb,rs
      REAL*8 cn,cnp,bnk,pb,rs
*
      REAL*8 fwhm,fwhm1,fwhm2,newt_raph,delta,f1,df1,f2,df2,fwhm_beam
      INTEGER i,itry
      LOGICAL warn
      EXTERNAL genfun,genfun2
*
      eps = 1.D-4
      looplim = -1
      ln2 = LOG(2.D0)
      pi = 4.D0*ATAN(1.D0)
      CALL initfac
* read FWHM for the beam
      READ(5,*) fwhm_beam
      pb = 4.D0*ln2/fwhm_beam**2
* read surface brightness profile
      np = 1
      rs = 1.D0
      DO i=1,PMAX
         READ(5,*,END=10) xp(np),yp(np)
         IF( xp(np).GE.0.D0 ) THEN
            IF( xp(np).GT.0.D0 .AND. yp(np).GT.0.D0 ) rs = xp(np)
            np = np+1
         ENDIF
      ENDDO
 10   np = np-1
* normalize to outer radius to prevent very large cn values
      DO i=1,np
         xp(i) = xp(i)/rs
      ENDDO
* get initial estimate for the FWHM...
      CALL getcn(1)
      itry = 1
 15   fwhm = 2.D0/(sqrt(cnp(0)/(cnp(1)*ln2)) - 0.36D0)
      IF( itry.EQ.2 ) fwhm = fwhm*1.5D0
      IF( itry.EQ.3 ) fwhm = fwhm/1.5D0
* and determine the final value
      delta = fwhm/10.D0
      i = 0
 20   CONTINUE
         IF( delta.LT.fwhm ) THEN
            fwhm1 = fwhm - delta
         ELSE
            fwhm1 = fwhm1/2.D0
         ENDIF
         call genfun(fwhm1,f1,df1,warn)
         fwhm2 = fwhm + delta
         call genfun(fwhm2,f2,df2,warn)
         IF( f1*f2.LT.0.D0 ) GOTO 30
         delta = 1.6D0*delta
         i = i+1
      IF( i.LE.20 ) GOTO 20
      WRITE(6,'("COULD NOT BRACKET ZERO")')
      STOP
 30   fwhm = newt_raph(genfun,fwhm1,f1,fwhm2,f2,eps,warn)
      itry = itry+1
      IF( warn .AND. itry.LE.3 ) GOTO 15
*
      IF( .NOT.warn ) WRITE(6,'("OBSERVED FWHM  = ",1P,E14.6)') fwhm*rs
*
      itry = 1
 35   fwhm = 2.D0/(sqrt(cn(0)/(cn(1)*ln2)) - 0.36D0)
      IF( itry.EQ.2 ) fwhm = fwhm*1.5D0
      IF( itry.EQ.3 ) fwhm = fwhm/1.5D0
* and determine the final value
      delta = fwhm/10.D0
      i = 0
 40   CONTINUE
         IF( delta.LT.fwhm ) THEN
            fwhm1 = fwhm - delta
         ELSE
            fwhm1 = fwhm1/2.D0
         ENDIF
         call genfun2(fwhm1,f1,df1,warn)
         fwhm2 = fwhm + delta
         call genfun2(fwhm2,f2,df2,warn)
         IF( f1*f2.LT.0.D0 ) GOTO 50
         delta = 1.6D0*delta
         i = i+1
      IF( i.LE.20 ) GOTO 40
      WRITE(6,'("COULD NOT BRACKET ZERO")')
      STOP
 50   fwhm = newt_raph(genfun2,fwhm1,f1,fwhm2,f2,eps,warn)
      itry = itry+1
      IF( warn .AND. itry.LE.3 ) GOTO 35
*
      IF( .NOT.warn ) WRITE(6,'("INTRINSIC FWHM = ",1P,E14.6)') fwhm*rs
*
      END
******************************************************************************
      SUBROUTINE initfac
*     
      INTEGER NLIM
      PARAMETER( NLIM=50 )
      COMMON/fac/ fac(0:NLIM)
      REAL*8 fac
      INTEGER i
      fac(0) = 1.D0
      DO i=1,NLIM
         fac(i) = DFLOAT(i)*fac(i-1)
      ENDDO
      RETURN
      END
******************************************************************************
      REAL*8 FUNCTION binom(n,k)
      INTEGER n,k
*     
      INTEGER NLIM
      PARAMETER( NLIM=50 )
      COMMON/fac/ fac(0:NLIM)
      REAL*8 fac
      binom = (fac(n)/fac(n-k))/fac(k)
      RETURN
      END
******************************************************************************
      SUBROUTINE getcn(i)
      INTEGER i
*     
      COMMON/parm/ ln2,pi,eps,looplim
      REAL*8 ln2,pi,eps
      INTEGER looplim
      INTEGER NLIM,PMAX
      PARAMETER( NLIM=50, PMAX=1024 )
      COMMON/fac/ fac(0:NLIM)
      REAL*8 fac
      COMMON/prof/ xp(PMAX),yp(PMAX),np
      REAL*8 xp,yp
      INTEGER np
      COMMON/cn/ cn(0:NLIM),cnp(0:NLIM),bnk(0:NLIM,0:NLIM),pb,rs
      REAL*8 cn,cnp,bnk,pb,rs
*
      INTEGER n,j,n22,n23
      REAL*8 c_norm,a,b,binom
      SAVE c_norm
*
      DO n=looplim+1,MIN(i,NLIM)
         cnp(n) = 0.D0
         n22 = 2*n + 2
         n23 = 2*n + 3
         DO j=2,np
            b = (yp(j)-yp(j-1))/(xp(j)-xp(j-1))
            a = yp(j) - b*xp(j)
            cnp(n) = cnp(n) + a/DFLOAT(n22)*(xp(j)**n22-xp(j-1)**n22) +
     $           b/DFLOAT(n23)*(xp(j)**n23-xp(j-1)**n23)
         ENDDO
         IF( n.EQ.0 ) c_norm = cnp(0)
         cnp(n) = cnp(n)/c_norm
         WRITE(6,'(I3,1P,E16.8)') n,cnp(n)
         cn(n) = cnp(n)
         DO j=0,n-1
            bnk(n,j) = binom(n,j)*(fac(n)/fac(j))
            cn(n) = cn(n) - bnk(n,j)/(pb*rs**2)**(n-j)*cn(j)
         ENDDO
*         WRITE(6,'(I3,1P,E16.8)') n,cn(n)
      ENDDO
      looplim = MIN(i,NLIM)
      RETURN
      END
******************************************************************************
      SUBROUTINE genfun(fwhm,genf,gendrv,warn)
      REAL*8 fwhm,genf,gendrv
      LOGICAL warn
*     
      INTEGER NLIM
      PARAMETER( NLIM=50 )
      COMMON/parm/ ln2,pi,eps,looplim
      REAL*8 ln2,pi,eps
      INTEGER looplim
      COMMON/fac/ fac(0:NLIM)
      REAL*8 fac
      COMMON/cn/ cn(0:NLIM),cnp(0:NLIM),bnk(0:NLIM,0:NLIM),pb,rs
      REAL*8 cn,cnp,bnk,pb,rs
      REAL*8 x,xt1,xt2,y
      INTEGER n,n2
*     
      warn = .FALSE.
      n = 0
      genf = 0.D0
      gendrv = 0.D0
 10   CONTINUE
         IF( n.GT.looplim ) CALL getcn(n)
         n2 = 2*n
         x = (-1.D0)**n
         y = 4.D0*ln2/fwhm**2
         xt1 = x*DFLOAT(n2+1)*cnp(n)*y**n/fac(n)
         IF( n.GT.0 ) THEN
            xt2 = -2.D0*x*DFLOAT(n2+1)*cnp(n)*y**n/fwhm/fac(n-1)
         ELSE
            xt2 = 0.D0
         ENDIF
         genf = genf + xt1
         gendrv = gendrv + xt2
         n = n+1
      IF( ABS(xt1).GT.ABS(gendrv*eps) .AND. n.LE.NLIM ) GOTO 10
*     
      IF( ABS(xt1).GT.ABS(gendrv*eps) ) THEN
         WRITE(0,'(''*** genfun: limit exceeded'')')
         warn = .TRUE.
      ENDIF
      IF( gendrv.GT.1.D10 .OR. gendrv.LT.0.D0 ) warn = .TRUE.
*     
      RETURN
      END
******************************************************************************
      SUBROUTINE genfun2(fwhm,genf,gendrv,warn)
      REAL*8 fwhm,genf,gendrv
      LOGICAL warn
*     
      INTEGER NLIM
      PARAMETER( NLIM=50 )
      COMMON/parm/ ln2,pi,eps,looplim
      REAL*8 ln2,pi,eps
      INTEGER looplim
      COMMON/fac/ fac(0:NLIM)
      REAL*8 fac
      COMMON/cn/ cn(0:NLIM),cnp(0:NLIM),bnk(0:NLIM,0:NLIM),pb,rs
      REAL*8 cn,cnp,bnk,pb,rs
      REAL*8 x,xt1,xt2,y
      INTEGER n,n2
*
      warn = .FALSE.
      n = 0
      genf = 0.D0
      gendrv = 0.D0
 10   CONTINUE
         IF( n.GT.looplim ) CALL getcn(n)
         n2 = 2*n
         x = (-1.D0)**n
         y = 4.D0*ln2/fwhm**2
         xt1 = x*DFLOAT(n2+1)*cn(n)*y**n/fac(n)
         IF( n.GT.0 ) THEN
            xt2 = -2.D0*x*DFLOAT(n2+1)*cn(n)*y**n/fwhm/fac(n-1)
         ELSE
            xt2 = 0.D0
         ENDIF
         genf = genf + xt1
         gendrv = gendrv + xt2
         n = n+1
      IF( ABS(xt1).GT.ABS(gendrv*eps) .AND. n.LE.NLIM ) GOTO 10
*     
      IF( ABS(xt1).GT.ABS(gendrv*eps) ) THEN
         WRITE(0,'(''*** genfun2: limit exceeded'')')
         warn = .TRUE.
      ENDIF
      IF( gendrv.GT.1.D10 .OR. gendrv.LT.0.D0 ) warn = .TRUE.
*     
      RETURN
      END
******************************************************************************
      REAL*8 FUNCTION newt_raph(func,x1,f1,x2,f2,xacc,warn)
      REAL*8 x1,f1,x2,f2,xacc
      LOGICAL warn
      EXTERNAL func
      INTEGER MAXIT
      PARAMETER( MAXIT=100 )
      INTEGER j
      REAL*8 df,dx,dxold,f,fh,fl,temp,xh,xl
      warn = .FALSE.
      fl = f1
      fh = f2
      IF( fl.EQ.0.D0 ) THEN
         newt_raph = x1
         RETURN
      ELSE IF( fh.EQ.0.D0 ) THEN
         newt_raph = x2
         RETURN
      ELSE IF( fl.LT.0.D0 ) then
         xl = x1
         xh = x2
      ELSE
         xh = x1
         xl = x2
      ENDIF
      newt_raph = 0.5D0*(x1+x2)
      dxold = ABS(x2-x1)
      dx = dxold
      CALL func(newt_raph,f,df,warn)
      DO j=1,maxit
         IF( ((newt_raph-xh)*df-f)*((newt_raph-xl)*df-f).GE.0.D0 .OR.
     *        ABS(2.D0*f).GT.ABS(dxold*df) ) THEN
            dxold = dx
            dx = 0.5D0*(xh-xl)
            newt_raph = xl+dx
            IF( xl.EQ.newt_raph ) RETURN
         ELSE
            dxold = dx
            dx = f/df
            temp = newt_raph
            newt_raph = newt_raph-dx
            IF( temp.EQ.newt_raph ) RETURN
         ENDIF
         IF( (xh-xl).LT.ABS(f/df) ) warn=.TRUE.
         IF( ABS(dx).LT.xacc ) RETURN
         CALL func(newt_raph,f,df,warn)
         IF( f.LT.0.D0 ) THEN
            xl = newt_raph
         ELSE
            xh = newt_raph
         ENDIF
      ENDDO
      warn = .TRUE.
      WRITE(0,'(''*** newt_raph: exceeding maximum iterations'')')
      RETURN
      END