C @(#)comcc.for	19.1 (ES0-DMD) 02/25/03 13:31:21
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
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENT:    COMCC			A. Lauberts     ESO	890711
C.KEYWORDS: table file, characteristic curve, image file
C.PURPOSE:  transform density D to intensity I according to the Llebaria formula
C           log I = C3*log(D-C1) + C4*log(C2-D) + C5
C           with calibration parameters C1 - C5
C.USAGE:    MIDAS applic program
C           execute as @@ COMCC P1 P2
C.PARAMETERS: P1  INA   C*8  table file with calibration parameters in 1st row
C             P2  OUTA	C*8  name of output image file with characteristic curve
C Output:   1-D image file with characteristic curve D = D(log I),
C           where log I limits are determined by descriptors NPIX,
C           START and STEP.
C.COMMENTS: derived from old routine CTAB
C---------------------------------------------------------------------------
	PROGRAM COMCC
C
	COMMON/VMR/MADRID(1)
	CHARACTER*60  INA,OUTA
	INTEGER      COL(9),STAT,PNTR
	REAL         CONST(9)
	LOGICAL      NULL(9)
C
	DATA COL/1,2,3,4,5,6,7,8,9/
	PARAMETER NCOL=9
C
C  set up MIDAS environment 
	CALL SXSPRO('COMCC',STAT)
C
C  get keywords
	CALL SXKGET('IN_A','C',1,60,IAV,INA,STAT)
	CALL SXKGET('OUT_A','C',1,60,IAV,OUTA,STAT)
C  read cc coeffs from table file
	CALL TXTOPN(INA,STAT)
	CALL TXRGET(INA,'R',1,NCOL,COL,CONST,NULL,STAT)
	IF(NULL(2).EQ..TRUE.) GOTO 999
	STEP=0.01
	START=-1.
	MAX=401
C  creat 1-D image
	CALL SXIPUT(OUTA,'O',1,MAX,START,STEP,OUTA,'LOG_INT',PNTR,NO,
     + STAT)
	CALL DENS(MADRID(PNTR),CONST,MAX,STEP,START)
C  write descriptor LHCUTS
	CALL SXDPUT(OUTA,'LHCUTS','R',CONST,1,4,STAT)
C  free data
999	CALL SXSEPI
	END
C
C
	SUBROUTINE DENS(A,CONST,MAX,STEP,START)
C
	REAL  A(MAX),CONST(9)
C
	F(X) = C3*ALOG10(X-C1) + C4*ALOG10(C2-X) + C5
	DF(X)= (C3/(X-C1) - C4/(C2-X)) / 2.3
C
	C1=CONST(2)
	C2=CONST(3)
	C3=CONST(4)
	C4=CONST(5)
	C5=CONST(6)
	D1=C1+.05
	D2=C2-.05
C  find D for log I = -1 to +3 using Newton-Raphson method
	D=D1
	FI=START
	FD=F(D)
	DFD=DF(D)
       DO I=1,MAX
	IT=0
10	D0=D
	IF(IT.GT.9) THEN
	  WRITE(6,*) ' >9 iterations at log I, D=',FI,D
	  RETURN
	ENDIF
	IT=IT+1
	F0=FD
	DF0=DFD
	D = D0 - (F0-FI) / DF0
	D = AMIN1(D2,AMAX1(D1,D))
	FD=F(D)
	DFD=DF(D)
	IF(ABS(D0-D).GT.0.0001.AND.ABS(FD-FI).GT.0.0001) GOTO 10
	A(I)=D
	FI=FI+STEP
       ENDDO
C  put LHCUTS in CONST
	CONST(1)=0.
	CONST(2)=D
	CONST(3)=0.
	CONST(4)=D
	RETURN
	END
C
