C @(#)dtom.for	17.1.1.1 (ES0-DMD) 01/25/02 17:19:20
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
	PROGRAM DTOM
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENT  program DTOM				A. Lauberts	ESO  890710
C
C.KEYWORDS  characteristic curve, Schmidt plates, standard galaxies,
C	    sensitometer spots
C
C.PURPOSE
C
C  Find least squares fit of
C  MPG = photographic magnitude/2.5 (in rings between circular apertures)
C                                    + innermost aperture
C  to
C  MEL = photoelectric magnitude/2.5 (in same rings as MPG)
C                                     curve fitted
C  and simultaneously
C
C   F  = A*log(D-DF) + B*log(DS-D) + CM (the D to log I function)
C  to  
C   E  = log I + CS (versus the spot densities D)
C
C  Minimize  S = Sum (MPG-MEL)**2 + Sum (F-E)**2
C      over   min0:max9 rings       min0:max21 spots
C  where 
C
C   D = Density measured by PDS (divide by 200 or 800!) 
C  DF = Density of the chemical fog  = X(1) = XX(1)  
C  DS = Density at saturation        = X(2) = XX(2)
C   G = Gamma = dD/d(log I)
C  RG = ln10/Gamma at infl           = X(3) = XX(7)       
C  DI = Density at inflexion point   = X(4) = XX(8)
C  DB = Density at sky background
C  CM = Zero point for mag/2.5       = X(5)
C       Zero point in sky units             = XX(5)
C  CS = Zero point for spots         = X(6) = XX(9)
C   A = RG*(DI-DF)**2/(DS-DF)               = XX(3)                     
C   B =-RG*(DS-DI)**2/(DS-DF)               = XX(4)
C  SM = Sky magnitude/20 micron aperture    = XX(6)
C
C.USAGE  MIDAS applic program
C	 execute as @@ DTOM
C
C.PARAMETERS
C
C Input: promted by inquire statements 
C
C  TABLEA    C*8	standard galaxy photoel table file
C			(format same as BVRLA.TBL used for the 
C			automated ESOLV procedures)
C  TABLEC    C*8	table file with calibration parameters
C  FRAMEA    C*8	image frame with central galaxy
C  FRAMEB    C*8	table file with sensitometer spots:
C  			SE   R*4 array  log I	intensity
C  			SD   R*4 array	D	density
C  NES	     R*4 array  1st elem, no of spots wanted
C  SPOTE     R*4	constant log I difference (option)
C  NA, NR    I*4	no of apertures (rings)
C  ESO	     R*4	ESO ident = (field.object no) = ESO.OBJ
C  SCAL	     R*4	plate scale arcseq/mm
C  WT	     R*4 array	weights of innermost, outermost, internal rings
C  SM	     R*4	sky mag/20 microns aperture (option)
C  XX	     R*4 array	initial calibration values (4 first elems)
C  FIX	     I*4 array	1/0 for fix/no-fix calibr parameters
C
C Descriptors:
C
C  BVR	     R*4 array 	ident ESO.OBJ + BVR vs log A coefficients read from 
C			TABLEA (format same as BVRLA.TBL used for the
C			automated ESOLV procedures).
C  GALCEN    R*4 array	Galaxy centre computed in program COMSKY
C			or set by interactive cursor
C  SKYBGR    R*4 array  Sky background computed in COMSKY
C			= command COMUTE/SKY
C
C Output:
C
C  CAL	     R*4 array	ident ESO.OBJ + calibration parameters -> TABLEC
C
C.COMMENTS
C
C Note: In case of non-existing V-R (all coeffs = 0's), a substitute
C	is taken as
C
C	V-R = 0.134 + 0.453 * (B-V)  at BS = 22 mags/sqarcs
C	d(V-R)/dlogA = -0.025 = constant,
C
C	being a good approximation (R.M.S. err = 0.02)
C	for 0.8 < B-V < 1.2 (E and S0's)
C 
C	The photoelectric values are reduced to outside the atmosphere
C	(airmass zero). Therefore, the present 'sky parameter' becomes 
C	brighter than the true sky. The magnitude difference is simply
C	the actual extinction. Whenever needed we take Bext = 0.20 (ESO
C	manual 1978) and Rext = 0.10 (ESO manual gives 0.08, own observ
C	0.12 after 1982) per unit airmass.
C
C	This program calls former NAGLIB subroutine E04GAF, 
C	which - together with supplementary smaller routines -
C	presently resides in module MNAG
C
C-------------------------------------------------------------------------------
	COMMON /VMR/MADRID(1)
	COMMON /DENS/ID(4096,9),MAXD,NR,WT(30),NS,DB,SM,A,B,XX,STEP,
     +  	SD(21),SE(21),EL(9),ORG,FIX
	DOUBLE PRECISION  RR,F,X(6),X02AAF,XTOL(6),SCALE(6),R(33),W(100),
     +		A,B,DB,ORG(6)
	REAL  XX(9),APER(9),MEL(9),RADA(9),WEIGHT(9),CEL(9),
     +		STEPA(2),STARTA(2),CEN(2),LOGA(9),BVR(31),CAL(9),SEL(9)
	INTEGER  NPIXA(2),STAT,PNTRA,NES(2),COLA(31),COLC(9),FIX(6)
	CHARACTER*1   CALC,CHR
	CHARACTER*60  FRAMEA,FRAMEB,SKYMAG,TABLEA,TABLEC,LABEL(2)
	CHARACTER*20  IDENTA
	CHARACTER*24  CUNITA
	CHARACTER*9   FSSQ
	CHARACTER*72  CALIB,RESIDS,SURFB
	CHARACTER*3   ERRNO,RIN,SPO
	LOGICAL	      NULL(31),IFL
C
	PARAMETER     NCOLA=31, NCOLC=9
	DATA          COLC/1,2,3,4,5,6,7,8,9/,
     +		      LABEL/'DENSITY','LOG_INT'/
C
	EXTERNAL  RESID,LSQ,MONIT
C
	G(D)=C3*ALOG10(D-C1)+C4*ALOG10(C2-D)+C5
C
	DO I=1,31
	  COLA(I)=I
	ENDDO
C
C  set up MIDAS environment, access tables (already existing)
C
	CALL SXSPRO('DTOM',STAT)
	CALL SXKGET('TABLEA','C',1,60,IAV,TABLEA,STAT)
	CALL SXKGET('TABLEC','C',1,60,IAV,TABLEC,STAT)
	CALL TXTOPN(TABLEA,STAT)
	CALL TXTOPN(TABLEC,STAT)
C  get plate info
	CALL SXKGET('INPUTC','C',1,1,IAV,CHR,STAT)
	CALL SXKGET('SCALE','R',1,1,IAV,SCAL,STAT)
	SKY=0.
	DB=0.
	DO 40 J=1,30
	WT(J)=1.
40	CONTINUE
	CALL SXKGET('IN_B','C',1,60,IAV,FRAMEB,STAT)
	CALL TXTOPN(FRAMEB,STAT)
	CALL TXIGET(FRAMEB,NCOLB,NROWB,IAV,STAT)
	NS=NROWB
	DO I=1,NS
	  CALL TXRGET(FRAMEB,'R',I,2,COLA,XX,NULL,STAT)
	  SE(I)=XX(1)
	  SD(I)=XX(2)
	ENDDO
	WRITE(6,111) (SD(I),I=1,NS)
111	FORMAT(' SPOTD',10F7.3)
48	WRITE(6,*) ' first element, no of spots: '
	READ(5,*) NES
	NE=NES(1)
	NS=NES(2)
	IF(NS.EQ.1.OR.NS.GT.21) THEN
	  CALL SXTPUT('NO OF SPOTS MUST BE 0 OR 2-21',STAT)
	  GOTO 48
	ENDIF
	IF(NS.EQ.0) GOTO 55
	WRITE(6,*) ' const log intensity difference 
     + (if =0 then 1st col in Cxxx): '
	READ(5,*) SPOTE
49	CONTINUE
	IF(NE.NE.1) THEN
	  NE=NE-1
	  DO I=1,NS
	    SD(I)=SD(I+NE)
	    SE(I)=SE(I+NE)
	  ENDDO
	ENDIF
	IF(SPOTE.GT.0.) THEN
	  DO I=1,NS
	    SE(I)=SPOTE*FLOAT(I)
	  ENDDO
	  S=-0.5
	  DO I=1,NROWB
	    CALL TXEPUT(FRAMEB,'R',I,1,S,STAT)
	    S=S+SPOTE
	  ENDDO
	ENDIF
	WRITE(6,112) (SE(I),I=1,NS)
112	FORMAT(' SPOTE',10F7.3)
	CALL TXTCLO(FRAMEB,STAT)
	NR=0
55	CALL SXKGET('NAPE','I',1,1,IAV,NA,STAT)
	IF(NA.GT.9) THEN
	  CALL SXTPUT(' >9 APERTURES/OBJECT',STAT)
	  GOTO 55
	ENDIF
C  get image frame 
	CALL SXKGET('IN_A','C',1,60,IAV,FRAMEA,STAT)
	CALL SXIGET(FRAMEA,'I',2,NAXISA,NPIXA,STARTA,STEPA,
     +              IDENTA,CUNITA,PNTRA,STAT)
	IF(NA.NE.0) CALL SXDGET(FRAMEA,'GALCEN','R',1,2,IAV,CEN,STAT)
	CALL SXDGET(FRAMEA,'SKYBGR','R',1,1,IAV,SKY,STAT)
	CALL SXDGET(FRAMEA,'O_TIME','R*8',1,1,IAV,A,STAT)
	IF(A.GT.1983.5D0) THEN
	  MAXD=4096
	  DB=SKY/800.
	 ELSE
	  MAXD=1024
	  DB=SKY/200.
	ENDIF
C  convert sky to mags/(20*20) square micron aperture
	STEP=STEPA(1)
	IF(NA.EQ.0) THEN
	  CALL SXKGET('SKYMAG','R',1,1,IAV,SM,STAT)
	  SM=SM-5.*ALOG10(SCAL*0.020)
	ENDIF
	NDIM=NPIXA(1)*NPIXA(2)
	IF(NA.EQ.0) GOTO 100
	NA1=9
	S=5000./SCAL
	A1=10.**0.1
	DO I=1,NA1
	  RADA(I)=S
	  S=S*A1
	ENDDO
	A1=0.9
	M=100*MAXD/1024
	CALL INTEGR(NDIM,MADRID(PNTRA),STARTA,STEPA,NPIXA,CEN,RADA,
     +              NA1)
	CALL SXTPUT(' D-UP    0.5  1.0  1.5  2.0  2.5  3.0  3.5  
     +4.0  4.5  5.0',STAT)
	DO J=1,NA1
	  DO I=1,M*10,M
	    DO K=1,M-1
	      ID(I,J)=ID(I,J)+ID(I+K,J)
	    ENDDO
	  ENDDO
	  A1=A1+0.1
	  WRITE(SURFB,204) A1, (ID(I,J), I=1,M*10,M)
	  CALL SXTPUT(SURFB,STAT)
	ENDDO
	S=RADA(NA1)
	NA1=5
	A1=10.**0.1
	DO I=1,NA1
	  RADA(I)=S
	  S=S*A1
	ENDDO
	A1=1.8
	CALL INTEGR(NDIM,MADRID(PNTRA),STARTA,STEPA,NPIXA,CEN,RADA,
     +              NA1)
	DO J=2,NA1
	  DO I=1,M*10,M
	    DO K=1,M-1
	      ID(I,J)=ID(I,J)+ID(I+K,J)
	    ENDDO
	  ENDDO
	  A1=A1+0.1
	  WRITE(SURFB,204) A1, (ID(I,J), I=1,M*10,M)
	  CALL SXTPUT(SURFB,STAT)
	ENDDO
	WRITE(6,*) ' first, last log aperture: '
	READ(5,*) LOGA(1),LOGA(2)
C  specify table entry
60	CALL SXKGET('ESOOBJ','R',1,1,IAV,ESO,STAT)
	APER(1)=10.**LOGA(1)
	RADA(1)=APER(1)*500./SCAL
       IF(NA.GT.1) THEN
	 S=(LOGA(2)-LOGA(1))/FLOAT(NA-1)
	 A1=10.**S
	DO I=2,NA
	 APER(I)=APER(I-1)*A1
	 LOGA(I)=LOGA(I-1)+S
	 RADA(I)=RADA(I-1)*A1
	ENDDO
       ENDIF
C  store density distribution in common block /DENS/ID
C  new RADA out!
	CALL INTEGR(NDIM,MADRID(PNTRA),STARTA,STEPA,NPIXA,CEN,RADA,
     +              NA)
C  refine apertures considering actual no of pixels
	DO I=1,NA
	  APER(I)=RADA(I)*SCAL*0.002
	  LOGA(I)=ALOG10(APER(I))
	ENDDO
C  read BVR vs log A polynomial fit coefficients 
C  BVR(1) = ESO ident in the fld.obj no form 
C  Example: 271.012 means ESO field = 271, 12 = object no in ESO/Upp cat
C
	ERR=0.0005
	CALL TXESER(TABLEA,'R*4',1,ESO,ERR,1,NSEQ,STAT)
	CALL TXRGET(TABLEA,'R*4',NSEQ,NCOLA,COLA,BVR,NULL,STAT)
C  compute array of U,B,V or R values vs aperture A
C  as well as surface brightness in mags/sqarcsec
	CALL POLY('B',BVR,LOGA,MEL,SEL,NA)
	CALL SXTPUT('B SURFBR MAGS/SQARCSEC',STAT)
	WRITE(SURFB,201) (SEL(I), I=1,NA)
	CALL SXTPUT(SURFB,STAT)
	IF(CHR.NE.'B') THEN
	  CALL POLY(CHR,BVR,LOGA,MEL,CEL,NA)
	  CALL SXTPUT(' B-'//CHR,STAT)
	  DO I=1,NA
	    CEL(I)=SEL(I)-CEL(I)
	  ENDDO
	  WRITE(SURFB,201) (CEL(I), I=1,NA)
	  CALL SXTPUT(SURFB,STAT)
	ENDIF
	CALL SXTPUT(' AT APERTURES ARCSEC',STAT)
	WRITE(SURFB,201) (APER(I), I=1,NA)
	CALL SXTPUT(SURFB,STAT)
C  take .4 times magnitude of inner aperture
	EL(1)=.4*MEL(1)
C  transform magnitudes to intensities
	DO 80 I=1,NA
80	MEL(I)=10.**(-.4*MEL(I))
	IF(NA.EQ.1) GOTO 100
C  take negative log of intensities in rings bound by apertures
	DO 90 I=2,NA
90	EL(I) = - ALOG10(MEL(I)-MEL(I-1))
100	NR=NA
C  free this area for next input
	CALL SXFCLO(FRAMEA,STAT)
120	N=6
	M=NR+NS
	IF(NR.EQ.0) GOTO 125
	CALL SXKGET('WEIGHT','R',1,9,IAV,WEIGHT,STAT)
	WT(1)=WEIGHT(1)
	WT(NA)=WEIGHT(2)
	DO I=2,NA-1
	  WT(I)=WEIGHT(I+1)
	ENDDO
C  equal importance photoel vs photographic
	S=1.
	IF(NS.NE.0) S=NS/FLOAT(NR)
	DO I=1,NR
	  WT(I)=WT(I)*S
	ENDDO
C  get initial values XX
125	CALL SXKGET('INPUTR','R',1,4,IAV,XX,STAT)
	CALL SXKGET('INPUTI','I',1,N,IAV,FIX,STAT)
      DO I=1,4
	X(I)=XX(I)
      ENDDO
	X(5)=0.
	X(6)=0.
	CALL RESID(M,N,X,R,IFL)
	IF(NR.EQ.0) GOTO 130
      A1=0.
      DO I=1,NR
	A1=A1+WT(I)
	X(5)=X(5)-R(I)
      ENDDO
	X(5)=X(5)/A1
130	CONTINUE
      A1=0.
      DO I=NR+1,M
	A1=A1+WT(I)
	X(6)=X(6)-R(I)
      ENDDO
	IF(NS.NE.0) X(6)=X(6)/A1
        XTOL(1)=DSQRT(X02AAF(RR))
      DO I=1,N
	ORG(I)=X(I)
      ENDDO
      DO I=2,N
	XTOL(I)=XTOL(1)
      ENDDO
        METHOD=1
        IW=N*(N+4)+M    
        IPRINT=1    
        MAXCAL=20      
        IFAIL=0
	CALL SXTPUT('E04GAF RESULTS:',STAT)
        CALL E04GAF(M,N,X,R,F,XTOL,METHOD,SCALE,W,IW,RESID,LSQ,   
     +              MONIT,IPRINT,MAXCAL,IFAIL)
160	WRITE(ERRNO,202) IFAIL
	CALL SXTPUT('THIS HAS ERROR NUMBER'//ERRNO,STAT)
C  update calibration parameters in TABLEC (one row only)
	CALL SXKGET('ESOFLD','I',1,1,IAV,NSEQ,STAT)
	IF(NA.NE.0.) THEN
	  CAL(1)=ESO
	  ELSE
	  CAL(1)=NSEQ
	ENDIF
	I=ESO
	IF(I.NE.NSEQ) TYPE *,ESO,' not in field ',NSEQ
       DO I=1,8
	CAL(I+1)=XX(I)
       ENDDO
	CALL TXRPUT(TABLEC,'R',1,NCOLC,COLC,CAL,STAT)
999	CALL TXTCLO(TABLEA,STAT)
	CALL TXTCLO(TABLEC,STAT)
	CALL SXSEPI
200	FORMAT(E9.3)
201	FORMAT(9F8.3)
202	FORMAT(I3)
203	FORMAT(F8.2)
204	FORMAT(1X,F3.1,3X,10I5)
	END                                                     
C   
C    
      SUBROUTINE RESID(M,N,XC,RC,IFL)     
C  CALLED BY EO4GAF AND LSQ 
C  CALCULATES THE VALUES OF THE RESIDUALS RC AT XC
     	COMMON /DENS/ID(4096,9),MAXD,NR,WT(30),NS,DB,SM,A,B,XX,STEP,
     +		SD(21),SE(21),EL(9),ORG,FIX
      LOGICAL IFL        
      DOUBLE PRECISION  XC,RC,Q,X41,X24,XC1,XC2,XC4,X,A,B,IP,DB,ORG(6)
      DIMENSION XC(N),RC(M),XX(9)
      INTEGER FIX(6)
C
	F(X) = A*DLOG10(X-XC1) + B*DLOG10(XC2-X)
	DO I=1,N
	  IF(FIX(I)) XC(I)=ORG(I)
	ENDDO
	XC1=XC(1)
	XC2=XC(2)
	XC4=XC(4)
	Q=XC(3)/(XC2-XC1)
	X41=XC4-XC1
	X24=XC2-XC4
      IF(XC1.GT.DB*.1) GOTO 3
      IF(X41.GT.0.D0.AND.X24.GT.0.D0) GOTO 4
	GOTO 3
4	A= Q*X41**2
	B=-Q*X24**2
	IF(NR.EQ.0) THEN
	  XC(5)=0.D0
	  GOTO 11
	ENDIF
      DO 1 I=1,NR
	CALL INTRING(XC,I,IP,IFL)
	IF(IFL.EQ..TRUE.) RETURN
1     RC(I) = ( DLOG10(IP) + EL(I) + XC(5) ) * WT(I)
11	CONTINUE
      DO 2 I=1,NS
	X=SD(I)
        IF(X.GT.XC1.AND.X.LT.XC2) GOTO 2
	GOTO 3
2     RC(I+NR) = ( F(X) - SE(I) +XC(6) ) * WT(I+NR)
        IFL=.FALSE.
        RETURN
3	IFL=.TRUE.
	CALL SXTPUT('IFL=.TRUE. FROM RESID',STAT)
	RETURN              
        END                                                                  
C  
C   
      SUBROUTINE LSQ(M,N,XC,RC,AJTJC,GC)    
C  CALCULATES THE JACOBIAN RC AT XC   
C  CALLED FROM E04GAF    
      DOUBLE PRECISION  XC,RC,AJTJC,GC,AJC,SUM,DELTA,DELTA2,XJ,YFIT
      DIMENSION XC(N),RC(M),GC(N),AJTJC(N,N),AJC(30,6),YFIT(30)     
      LOGICAL IFL
      PARAMETER (DELTA=.001D0)
      PARAMETER (DELTA2=.002D0)
C  NON ANALYTICAL DERIVATIVE
      DO 3 J=1,N
      XJ=XC(J)
      XC(J)=XJ+DELTA
      CALL RESID(M,N,XC,RC,IFL)
      DO 1 I=1,M
1     YFIT(I)=RC(I)
      XC(J)=XJ-DELTA
      CALL RESID(M,N,XC,RC,IFL)
      DO 2 I=1,M
2     AJC(I,J)=(YFIT(I)-RC(I))/DELTA2
3     XC(J)=XJ
      DO 80 I=1,N    
         SUM=0.    
         DO 20 K=1,M   
            SUM=SUM+AJC(K,I)*RC(K)     
20       CONTINUE     
         GC(I)=SUM      
         II=I    
         DO 60 J=1,II   
            SUM=0.  
            DO 40 K=1,M  
               SUM=SUM+AJC(K,I)*AJC(K,J)  
40          CONTINUE  
            AJTJC(J,I)=SUM  
60       CONTINUE  
80    CONTINUE     
      RETURN   
      END    
C     
C   
      SUBROUTINE MONIT(M,N,XC,RC,FC,GC,NCALL)     
C  CALLED BY E04GAF   
C  PRINTS THE VALUES EVERY PRINT ITERATIONS
	COMMON /DENS/ID(4096,9),MAXD,NR,WT(30),NS,DB,SM,A,B,XX,STEP
	DOUBLE PRECISION  XC,GC,RC,FC,DB,A,B                         
	DIMENSION XC(N),GC(N),RC(M),R(30),XX(9)
	CHARACTER*3   CALL,RIN,SPO
	CHARACTER*9   SSQ
	CHARACTER*72  CALIB,RESIDS
	WRITE(CALL,200) NCALL
200	FORMAT(I3)
	WRITE(SSQ,201) FC
201	FORMAT(E9.3)
	CALL SXTPUT('AFTER'//CALL//' CALLS OF RESID, THE SUM OF
     + SQUARES IS '//SSQ//' AT',STAT)
	CALL SXTPUT('   DF      DS      A       B       CM
     +      SM      RG      DI',STAT)
      DO 1 I=1,5
1     XX(I)=XC(I)
        XX(3)=A
        XX(4)=B
	XX(9)=XC(6)
C  scale intensities to sky units through adjusting XX(5)
	S=DB
	IF(S.LE.XX(1)) THEN
	  WRITE(6,*) ' SKYBGR < FOG!'
	  RETURN
	ENDIF
	S=A*ALOG10(S-XX(1))+B*ALOG10(XX(2)-S)+XX(5)
	XX(5)=XX(5)-S
C  sky magnitude/20 micron aperture
	IF(NR.NE.0) THEN
	  XX(6)=-2.5*(S - 2.*ALOG10(STEP/20.))
	 ELSE
	  XX(6)=SM
	ENDIF
	XX(7)=XC(3)
	XX(8)=XC(4)
	WRITE(CALIB,202) (XX(I),I=1,8)
202	FORMAT(9F8.3)
	CALL SXTPUT(CALIB,STAT)
	WRITE(RIN,200) NR
	WRITE(SPO,200) NS
	CALL SXTPUT('WITH RESIDUES OF'//RIN//' RINGS AND'
     +//SPO//' SPOTS',STAT)
      DO I=1,M
	R(I)=RC(I)/WT(I)
      ENDDO
	IF(NR.NE.0) THEN
	  WRITE(RESIDS,202) (R(I), I=1,NR)
	  CALL SXTPUT(RESIDS,STAT)
	ENDIF
	DO 20 I1=NR+1,37,9
	IF(M.LT.I1) GOTO 30
	I2=I1+8
	IF(M.LT.I2) I2=M
	WRITE(RESIDS,202) (R(I), I=I1,I2)
20	CALL SXTPUT(RESIDS,STAT)
30	CONTINUE
      RETURN
      END
C
C
	SUBROUTINE INTRING(X,IR,IP,IFL)
C  called from RESID
C  computes integrated intensities in rings
	COMMON /DENS/ID(4096,9),MAXD,NR,WT(30),NS,DB,SM,A,B
	DOUBLE PRECISION  X(6),A,B,DB,DF,DF1,DS,DS1,DK,DI,IP,IB,IT
	LOGICAL IFL
C
	IF(MAXD.EQ.4096) THEN
	  DI=0.00125D0
	ELSE
	  DI=0.005D0
	ENDIF
	DF=X(1)
	DS=X(2)
	DF1=DF+.01D0
	DS1=DS-.01D0
	IF(DB.GT.DF1.AND.DB.LT.DS1) GOTO 4
	GOTO 2
4       IB = (DB-DF)**A * (DS-DB)**B
	IP=0.
	DK=0.
      DO 1 I=1,MAXD
	II=ID(I,IR)
	IF(II.EQ.0) GOTO 1
	IF(DK.GT.DF1.AND.DK.LT.DS1) GOTO 3
	GOTO 2
3	IT = (DK-DF)**A * (DS-DK)**B
        IP = IP + (IT-IB)*FLOATJ(II)
1     DK=DK+DI
	IFL=.FALSE.
	IF(IP.GT.0.D0) RETURN
2	IFL=.TRUE.
	CALL SXTPUT('IFL=.TRUE. FROM INTRING',STAT)
	RETURN
	END
C
C
	SUBROUTINE INTEGR(NDIM,AA,STARTA,STEPA,NPIXA,CEN,RADA,NA)
C  called from DTOM
C  calculates distribution of densities in ring-shaped regions
	COMMON /DENS/ID(4096,9),MAXD
	REAL  AA(NDIM),STARTA(2),STEPA(2),CEN(2),RADA(9),RAD2(9)
	INTEGER   NPIXA(2),NUMR(9)
	DO J=1,NA
	  DO I=1,MAXD
	    ID(I,J)=0
	  ENDDO
	  NUMR(J)=0
	ENDDO
      DO I=1,NA
	RAD2(I)=RADA(I)**2
      ENDDO
C  put restrictions on max region to search
	X0=CEN(1)
	Y0=CEN(2)
	RAD=RADA(NA)
	X1=X0-RAD
	X2=X0+RAD
	Y1=Y0-RAD
	Y2=Y0+RAD
	S=1.
	J1=(Y1-STARTA(2))/STEPA(2) -S
	J2=(Y2-STARTA(2))/STEPA(2) +S +2.
	IF(J1.LT.1) J1=1
	IF(J2.GT.NPIXA(2)) J2=NPIXA(2)
	IF(J1.GE.J2) RETURN
	Y1=STARTA(2)+(J1-1)*STEPA(2)
	I1=(X1-STARTA(1))/STEPA(1) -S
	I2=(X2-STARTA(1))/STEPA(1) +S +2.
	IF(I1.LT.1) I1 = 1
	IF(I2.GT.NPIXA(1)) I2 = NPIXA(1)
	IF(I1.GE.I2) RETURN
	X1=STARTA(1)+(I1-1)*STEPA(1)
	Y=Y1
	DO 3 J = J1,J2
	LOC=(J-1)*NPIXA(1)+I1
	YY=(Y-Y0)**2
	X=X1
	DO 2 I = I1,I2
	D2=(X-X0)**2 + YY
	DO 4 IR=1,NA
	IF(D2.GT.RAD2(IR)) GOTO 4
	IL=JIFIX(AA(LOC)+1.)
	IL=JMIN0(IL,MAXD)
C  form distribution of densities in ringshaped regions
	ID(IL,IR)=ID(IL,IR)+1
	NUMR(IR)=NUMR(IR)+1
	GOTO 1
4	CONTINUE
1	LOC=LOC+1
2	X=X+STEPA(1)
3	Y=Y+STEPA(2)
C  refine radii considering actual no of pixels
	IF(NA.GT.1) THEN
	  DO I=NA,2,-1
	    DO J=1,I-1
		NUMR(I)=NUMR(I)+NUMR(J)
	    ENDDO
	  ENDDO
	ENDIF
	DO I=1,NA
	  RADA(I)=SQRT(NUMR(I)/3.141592)*STEPA(1)
	ENDDO
	RETURN
	END
C
	SUBROUTINE POLY ( T,C,LOGA,MEL,SEL,NA )
C  compute magnitudes MEL and mean surface brightess SEL
C
	CHARACTER*1  T
	REAL	     C(1),MEL(NA),SEL(NA),LOGA(9)
	DO I=1,NA
	  X=LOGA(I)
	  IF(T.EQ.'B') THEN
	    IF(C(18).EQ.0.) THEN
		C(14)=-0.3456+C(20)*2.2502
		C(15)=-0.0344+C(21)*2.1358
		C(16)= 0.0451+C(22)*2.1070
		C(17)=-0.0143+C(23)*2.0661
	    ENDIF
	    MEL(I)=C(2)+C(14) + X*(C(3)+C(15) + X*(C(4)+C(16) +
     +		   X*(C(5)+C(17))))
	   ELSE IF(T.EQ.'R') THEN
	    IF(C(24).EQ.0.) THEN
		C(20)= 0.1536+C(14)*0.4444
		C(21)= 0.0161+C(15)*0.4682
		C(22)=-0.0214+C(16)*0.4746
		C(23)= 0.0069+C(17)*0.4840
	    ENDIF
	    MEL(I)=C(2)-C(20) + X*(C(3)-C(21) + X*(C(4)-C(22) +
     +		   X*(C(5)-C(23))))
	   ELSE IF(T.EQ.'U') THEN
	    MEL(I)=C(2)+C(8)+C(14)+X*(C(3)+C(9)+C(15)+
     +             X*(C(4)+C(10)+C(16)+X*(C(5)+C(11)+C(17))))
	   ELSE IF(T.EQ.'V') THEN
	    MEL(I)=C(2)+X*(C(3)+X*(C(4)+X*C(5)))
	   ELSE
	    CALL SXTPUT(' LEADING CHAR MUST BE U, B, V or R !',STAT)
	    RETURN
	  ENDIF
	  IF(I.NE.1.AND.MEL(I).LT.MEL(I1)) THEN
	    SEL(I)= MEL(I) +5.*X -0.262
     +		   -2.5*ALOG10(1.-10.**(-.4*(MEL(I1)-MEL(I))))
     +		   +2.5*ALOG10(1.-10.**(2.*(X1-X)))
	   ELSE
	    SEL(I)=MEL(I) +5.*X -0.262
	  ENDIF
	  I1=I
	  X1=X
	ENDDO
	END