program combine
!----------------------------------------------------------------------
! GILDAS Combine in different ways two input images
! (or data cubes)...
!----------------------------------------------------------------------
use image_def
use gkernel_interfaces
use gbl_format
!
character(len=filename_length) :: 80 namex,namey,namez
character(len=20) :: code
logical error
real ay,az,ty,tz,b,c
type (gildas) :: hx, hy, hz
real, allocatable :: dx(:,:), dy(:,:), dz(:)
integer(kind=4) :: i, j, n, m
integer :: ier
!
call gildas_open
call gildas_char('Z_NAME$',namez)
call gildas_real('Z_FACTOR$',az,1)
call gildas_real('Z_MIN$',tz,1)
call gildas_char('Y_NAME$',namey)
call gildas_real('Y_FACTOR$',ay,1)
call gildas_real('Y_MIN$',ty,1)
call gildas_char('X_NAME$',namex)
call gildas_real('BLANKING$',b,1)
call gildas_real('OFFSET$',c,1)
call gildas_char('FUNCTION$',code)
call gildas_close
!
n = len_trim(namez)
if (n.eq.0) goto 100
call gildas_null(hz)
call sic_parsef(namez(1:n),hz%file,' ','.gdf')
call gdf_read_header(hz,error)
if (error) then
call gagout('F-COMBINE, Cannot read input file')
goto 100
endif
!
n = len_trim(namey)
if (n.eq.0) goto 100
call gildas_null(hy)
call sic_parsef(namey(1:n),hy%file,' ','.gdf')
call gdf_read_header(hy,error)
if (error) then
call gagout('F-COMBINE, Cannot read input file')
goto 100
endif
!
if (hz%gil%eval.ge.0.0) hz%gil%eval = &
max(hz%gil%eval,abs(hz%gil%bval*1e-7))
if (hy%gil%eval.ge.0.0) hy%gil%eval = &
max(hy%gil%eval,abs(hy%gil%bval*1e-7))
!
! Check input dimensions
do i=1,gdf_maxdims
if (hy%gil%dim(i).ne.hz%gil%dim(i)) then
n = 1
do j=i,gdf_maxdims
n = n*hz%gil%dim(j)
enddo
if (n.ne.1) then
call gagout('F-COMBINE, Input images are non coincident')
goto 100
else
call gagout('W-COMBINE, Combining a cube with a plane')
endif
endif
enddo
!
call gdf_copy_header(hy,hx)
n = len_trim(namex)
if (n.eq.0) goto 100
call sic_parsef(namex(1:n),hx%file,' ','.gdf')
hx%gil%blan_words = 2
hx%gil%bval = b
hx%gil%eval = 0.0
hx%gil%extr_words = 0 ! No extrema computed
!
! Allocate the arrays. Note that the allocated arrays do not conform
! to the shape of the images: DZ is allocated as a 1-D array, DX,DY
! as 2-D arrays, possibly of second dimension 1.
!
n = hz%loca%size
m = hx%loca%size/hz%loca%size
allocate(dx(n,m),dy(n,m),dz(n),stat=ier)
if (ier.ne.0) then
call gagout('F-COMBINE, Input images are non coincident')
goto 100
endif
!
! Read the input data
call gdf_read_data(hz,dz,error)
call gdf_read_data(hy,dy,error)
!
if (code.eq.'ADD') then
call add002(dz,dy,dx, &
n,m, &
hz%gil%bval,hz%gil%eval,az,tz, &
hy%gil%bval,hy%gil%eval,ay,ty, &
hx%gil%bval,c)
elseif (code.eq.'DIVIDE') then
call div002(dz,dy,dx, &
n,m, &
hz%gil%bval,hz%gil%eval,az,tz, &
hy%gil%bval,hy%gil%eval,ay,ty, &
hx%gil%bval,c)
elseif (code.eq.'MULTIPLY') then
call mul002(dz,dy,dx, &
n,m, &
hz%gil%bval,hz%gil%eval,az,tz, &
hy%gil%bval,hy%gil%eval,ay,ty, &
hx%gil%bval,c)
elseif (code.eq.'OPTICAL_DEPTH') then
call opt002(dz,dy,dx, &
n,m, &
hz%gil%bval,hz%gil%eval,az,tz, &
hy%gil%bval,hy%gil%eval,ay,ty, &
hx%gil%bval,c)
else
call gagout('Invalid operation code '//code)
goto 100
endif
!
! Write ouput file
call gdf_write_image(hx,dx,error)
!
stop 'S-COMBINE, Successful completion'
!
100 call sysexi (fatale)
end
!
subroutine add002(z,y,x,n,m,bz,ez,az,tz,by,ey,ay,ty,bx,c)
!---------------------------------------------------------------------
! GDF Internal routine
! Linear combination of input arrays
! X = Ay*Y + Az*Z + C
! Arguments
! Z R*4(*) Input array (N)
! Y R*4(*) Input array (N,M)
! X R*4(*) Output array (N,M)
! N,M I Dimensions of arrays
! BX,BY,BZ R*4 Blanking values
! EY,EZ R*4 Tolerance on blanking
! AZ,AY R*4 Multiplicative factor of array Z, Y
! TZ,TY R*4 Threshold on Z,Y
! C R*4 Additive constant
!---------------------------------------------------------------------
integer :: n !
real :: z(n) !
integer :: m !
real :: y(n,m) !
real :: x(n,m) !
real :: bz !
real :: ez !
real :: az !
real :: tz !
real :: by !
real :: ey !
real :: ay !
real :: ty !
real :: bx !
real :: c !
! Local
integer :: i,k
!
do k=1,m
do i=1,n
if (abs(z(i)-bz).gt.ez .and. abs(y(i,k)-by).gt.ey &
& .and. z(i).gt.tz .and. y(i,k).gt.ty) then
x(i,k) = ay*y(i,k) + az*z(i) + c
else
x(i,k) = bx
endif
enddo
enddo
end