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