C @(#)dtoi.for	19.1 (ES0-DMD) 02/25/03 13:31:22
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 DTOI
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENT
C
C  program  DTOI			A. Lauberts     ESO	890711
C
C.KEYWORDS 
C
C  Density, intensity, calibrated Schmidt plates
C
C.PURPOSE
C
C  transform density D to intensity I according to the Llebaria formula
C
C	log I = C3*log(D-C1) + C4*log(C2-D) + C5
C
C  with calibration parameters C1 - C5 read from a one row table file
C
C.USAGE  MIDAS applic program
C	 execute as @@ DTOI P1 P2
C
C.PARAMETERS
C
C P1	FRAMEA	C*8	name of input image frame in PDS density units
C P2	TABLEC	C*8	name of table file with calibration parameters
C P3	FRAMEC	C*8	name of output image frame in calibrated intensity
C			scaled with sky background = 1
C
C.DESCRIPTORS (special attention)
C
C 	HISTORY C       checked for and appended with 'INTE' 
C 	O_TIME  R*8     year,month,day of PDS scan 
C			Note: Different prestanda before/after 1984.1.14
C			NOT automatically set by PDS; updated in COMSKY
C 	SKYBGR  R*4	array, computed in program COMSKY (mean and x,y terms,
C                   	4th element contains sky magnitude/20 micron aperture
C
C.RESTRICTIONS  Input digitized image in PDS units 
C		(new dynyamic range since 19841.14!)
C		Prepared by first running program COMSKY
C
C.COMMENTS  note that on output the 4th element in desc SKYBGR now contains
C	    sky mag/pixel !
C
C---------------------------------------------------------------------------
	IMPLICIT  INTEGER  (A-Z)
C
	COMMON/VMR/MADRID(1)
	CHARACTER*60     FRAMEA,FRAMEC,TABLEC
	CHARACTER*200    HISTORY
	CHARACTER*21     IDENTA
	CHARACTER*24     CUNITA
	INTEGER          NPIXA(2),COL(9),PNTRA,PNTRC
	REAL             STEPA(2),STARTA(2),SKY(4),CONST(9),ESO
	DOUBLE PRECISION TIME
	LOGICAL          NULL(9)
C
	DATA COL/1,2,3,4,5,6,7,8,9/
	PARAMETER NCOL=9
C
C  set up MIDAS environment 
C
	CALL SXSPRO('DTOI',STAT)
C
C  get keywords
C
	CALL SXKGET('IN_A','C',1,60,IAV,FRAMEA,STAT)
	CALL SXKGET('IN_B','C',1,60,IAV,TABLEC,STAT)
	CALL SXKGET('OUT_A','C',1,60,IAV,FRAMEC,STAT)
C
C  activate I/O images; note FRAMEC initially copy of FRAMEA
C
	CALL SXIGET(FRAMEA,'I',2,NAXISA,NPIXA,STARTA,STEPA,
     +              IDENTA,CUNITA,PNTRA,STAT)
	CALL SXIGET(FRAMEC,'IO',2,NAXISA,NPIXA,STARTA,STEPA,
     +              IDENTA,CUNITA,PNTRC,STAT)
C
C  get calibration parameters from table file, first row
C
	CALL TXTOPN(TABLEC,STAT)
	CALL TXRGET(TABLEC,'R',1,NCOL,COL,CONST,NULL,STAT)
C
C  read descriptors 
C
	CALL SXDGET(FRAMEA,'HISTORY','C',1,200,IAV,HISTORY,STAT)
	IF(INDEX(HISTORY,'DENSITY TO INTENSITY').NE.0) THEN
	  CALL SXTPUT(' ALREADY INTENSITY',STAT)
	  GOTO 999
	ENDIF
	CALL SXDGET(FRAMEA,'O_TIME','R*8',1,1,IAV,TIME,STAT)
	IF(TIME.GT.1983.5D0) THEN
	  MAXD=4096
	  SKY(1)=800.
	 ELSE
	  MAXD=1024
	  SKY(1)=200.
	ENDIF
	DO I=2,4
	  SKY(I)=0.
	ENDDO
	CALL SXDGET(FRAMEA,'SKYBGR','R',1,3,IAV,SKY,STAT)
C
	SIZE=NPIXA(1)*NPIXA(2)
C
C  do the D to I transformation
C
	CALL DENS(SIZE,MADRID(PNTRA),MADRID(PNTRC),CONST,SKY,MAXD)
C
	SKY(4)=SKY(4)-5.*ALOG10(STEPA(1)/20.)
	CALL SXDPUT(FRAMEC,'SKYBGR','R',SKY,1,4,STAT)
C  free data
999	CALL SXSEPI
	END
C
C
	SUBROUTINE DENS(NDIM,A,C,CONST,SKY,MAXD)
C
	REAL A(NDIM),C(NDIM),T(4096),CONST(9),SKY(4)
C
	F(X) = C3*ALOG10(X-C1) + C4*ALOG10(C2-X) + C5
C
	C1=CONST(2)
	C2=CONST(3)
	C3=CONST(4)
	C4=CONST(5)
	C5=CONST(6)
	D1=C1+.05
	D2=C2-.05
	IF(MAXD.EQ.4096) THEN
	  DI=0.00125
	 ELSE
	  DI=0.005
	ENDIF
	DO I=1,3
	  SKY(I)=SKY(I)*DI
	ENDDO
	FS=F(SKY(1))
C+
C  make look-up table to speed up transformation,
C  intensity in sky units
C-
       DO I=1,MAXD
	D = (I-1)*DI
	D = AMIN1(D2,AMAX1(D1,D))
	T(I) = 10.**(F(D)-FS)
       ENDDO
C
C  transform D to I
C
       DO N=1,NDIM
	C(N)=T(JINT(A(N)+1.))
       ENDDO
C
C  DIDD=dI/dD normalized to unit sky
C
	DIDD=(C3/(SKY(1)-C1)-C4/(C2-SKY(1)))/(10.**FS*ALOG(10.))
	SKY(1)=1.
	SKY(2)=SKY(2)*DIDD
	SKY(3)=SKY(3)*DIDD
	SKY(4)=CONST(7)-2.5*FS
	END
C
