C @(#)lsfun2i.for	19.1 (ES0-DMD) 02/25/03 13:17:40
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 @(#)lsfun2i.for	2.1.1.1 (ESO-IPG) 9/27/91 18:10:06
      SUBROUTINE LSFU2I(NRRES,NRPRM,PRM,FCTVEC,FCTJAC,LJAC,PIXEL,             
     +                    WEIGHT,IDIM1,IDIM2,IDIM3)                             
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++        
C                                                                               
C.MODULE                                                                        
C       FIT                                                                     
C                                                                               
C.NAME                                                                          
C       LSFUN2_I                                                                
C                                                                               
C.PURPOSE                                                                       
C   Function to compute the approximating function and the derivatives          
C   at each point for FIT/TABLE                                                 
C                                                                               
C.KEYWORDS                                                                      
C       Approximating Function.                                                 
C                                                                               
C.DESCRIPTION                                                                   
C       Test the type of FIT and go to respective routines                      
C                                                                               
C.LANGUAGE                                                                      
C       FORTRAN                                                                 
C                                                                               
C.CALLING SEQUENCE                                                              
C       CALL LSFUN2_I(NRRES,NRPRM,PRM,FCTVEC,FCTJAC,LJAC,                       
C    +                     PIXEL,WEIGHT,IDIM1,IDIM2,IDIM3)                      
C                                                                               
C.INPUT PARAMETERS                                                              
C       NRRES               INTEGER  Number of residuals                        
C       NRPRM             INTEGER  Number of parameters                         
C       PRM (NRPRM)       DOUBLE   Parameters                                   
C       LJAC                INTEGER  1st dimension of jacobian                  
C       PIXEL (IDIM_)       REAL     pixel values of the image                  
C       WEIGHT (IDIM_)      REAL     weighting factors                          
C       IDIM_               INTEGER  dimension of the 3-dim image.              
C                                                                               
C.MODIFIED PARAMETERS                                                           
C       none                                                                    
C                                                                               
C.OUTPUT PARAMETERS                                                             
C       FCTVEC (NRRES)      DOUBLE   Approximating values                       
C       FCTJAC (LJAC,NRRES) DOUBLE   values of the jacobian                     
C                                                                               
C.FILES                                                                         
C       FIT_NAG.INC/NOLIST                                                      
C                                                                               
C.MODULES CALLED                                                                
C       FTFUNC                                                                  
C                                                                               
C.AUTHOR                                                                        
C       Ph. DEFERT,      Feb 1986                                               
C                                                                               
C.MODIFICATIONS                                                                 
C        M Peron Oct 89 : bug fixed in the calculation of the derivatives 
C                Apr 92 : bug fixed in the calcul of coordinates
C                                                                               
C-----------------------------------------------------------------------        
C      IMPLICIT NONE                                                            
C     ..                                                                        
C     .. Scalar Arguments ..                                                    
      INTEGER LJAC,NRPRM,NRRES,IDIM1,IDIM2,IDIM3,MADRID(1)                      
C     ..                                                                        
C     .. Array Arguments ..                                                     
      DOUBLE PRECISION FCTJAC(LJAC,NRPRM),FCTVEC(NRRES),PRM(NRPRM)              
      REAL PIXEL(IDIM1,IDIM2,IDIM3),WEIGHT(IDIM1,IDIM2,IDIM3)                   
C     ..                                                                        
C     .. Scalars in Common ..                                                   
      INTEGER NRCOL,ISTAR                                                       
      CHARACTER WGTTYP*1                                                        
C     ..                                                                        
C     .. Arrays in Common ..                                                    
      INTEGER ICOL(10)                                                          
C     ..                                                                        
C     .. Local Scalars ..                                                       
      DOUBLE PRECISION W,Y,Y1,YOUT,YY                                           
      INTEGER IFUN,IPNT,K,BEGFCT,IP,I1,I2,I3,NFJAC                              
C     ..                                                                        
C     .. Local Arrays ..                                                        
      DOUBLE PRECISION DUMMY(10)                                                
      REAL X(10)                                                                
C     ..                                                                        
C     .. Common blocks ..                                                       
C     ..                                                                        
C     .. External Files ..                                                      
       INCLUDE 'MID_INCLUDE:FITNAGI.INC/NOLIST'                                 
      COMMON /LSQFUN/ICOL,NRCOL,ISTAR,WGTTYP                                    
      COMMON /VMR/MADRID                                                     
       INCLUDE 'MID_INCLUDE:FITNAGC.INC/NOLIST'               
C     ..                                                                        
C     .. Executable Statements ..                                               
C                                                                               
C   Deal with simple linear constraints                                         
C                                                                               
      DO 10 K = 1,NRPRM                                                         
          IP     = FIXPAR(K)                                                    
          IF (IP.EQ.0) PRM(K) = PARINI(K)                                       
          IF (IP.GT.0) THEN                                                     
              PRM(K) = PRM(IP)*PRPFAC(K)                                        
              PARAM(K) = PRM(K)                                                 
          END IF                                                                
                                                                                
   10 CONTINUE                                                                  
C                                                                               
C   Go through the image                                                        
C                                                                               
      X(1)   = START(1)                                                         
      X(2)   = START(2)                                                         
      X(3)   = START(3)                                                         
      IPNT   = 0                                                                
      DO 60 I3 = 1,NRPIX(3)                                                     
          X(2) = START(2)
          DO 50 I2 = 1,NRPIX(2)                                                 
              X(1) = START(1)
              DO 40 I1 = 1,NRPIX(1)                                             
C                                                                               
C  compute weights                                                              
C                                                                               
                  IPNT   = IPNT + 1                                             
                  Y      = PIXEL(I1,I2,I3)                                      
                  IF (WGTTYP(1:1).EQ.'C') THEN                                  
                      W      = 1.                                               
                                                                                
                  ELSE IF (WGTTYP(1:1).EQ.'W') THEN                             
                      W      = WEIGHT(I1,I2,I3)                                 
                                                                                
                  ELSE IF (WGTTYP(1:1).EQ.'S') THEN                             
                      YY     = ABS(Y)                                           
                      IF (YY.LT.1.E-12) THEN                                    
                          W      = 1.                                           
                                                                                
                      ELSE                                                      
                          W      = 1./YY                                        
                      END IF                                                    
                                                                                
                  ELSE IF (WGTTYP(1:1).EQ.'I') THEN                             
                      W      = 1./WEIGHT(I1,I2,I3)**2                           
                  END IF                                                        
                                                                                
                  Y1     = 0.D0                                                 
                  BEGFCT = 1                                                    
                  DO 30 IFUN = 1,NRFUN                                          
                      CALL FTFUNC(FCTCOD(IFUN),NRIND,X,ACTPAR(IFUN),            
     +                            PRM(BEGFCT),YOUT,DUMMY)                       
                      Y1     = Y1 + YOUT                                        
                      DO 20 NFJAC = 1,ACTPAR(IFUN)                              
                          K      = BEGFCT + NFJAC - 1                           
                          FCTJAC(IPNT,K) = SQRT(W)*DUMMY(NFJAC)                 
   20                 CONTINUE                                                  
                      BEGFCT = BEGFCT + ACTPAR(IFUN)                            
   30             CONTINUE                                                      
                  FCTVEC(IPNT) = SQRT(W)* (Y1-Y)                                
                  X(1)   = X(1) + STEP(1)                                       
   40         CONTINUE                                                          
              X(2)   = X(2) + STEP(2)                                           
   50     CONTINUE                                                              
          X(3)   = X(3) + STEP(3)                                               
   60 CONTINUE                                                                  
C                                                                               
C   Correct derivatives if linear constrained                                   
C                                                                               
      DO 90 K = 1,NRPRM                                                         
          IP     = FIXPAR(K)                                                    
          IF (IP.EQ.0) THEN                                                     
              DO 70 IPNT = 1,NRRES                                              
                  FCTJAC(IPNT,K) = 0.D0                                         
   70         CONTINUE                                                          
                                                                                
          ELSE IF (IP.GT.0) THEN                                                
              DO 80 IPNT = 1,NRRES                                              
                  IF (FIXPAR(IP).NE.0) THEN
                     FCTJAC(IPNT,IP) = FCTJAC(IPNT,IP) +                     
     +                              FCTJAC(IPNT,K)/PRPFAC(K)                    
                  ENDIF
                  FCTJAC(IPNT,K) = 0.D0                                     
   80         CONTINUE                                                          
          END IF                                                                
                                                                                
   90 CONTINUE                                                                  
      RETURN                                                                    
                                                                                
      END                                                                       

