next up previous contents index
Next: Initialization file Up: A Template Task Previous: A Template Task   Contents   Index

Source code

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



Gildas manager 2024-03-28