; $Id: tassb.pro,v 1.2 2002/01/14 20:59:27 amunter Exp $
function combine,a,ptnum,b,npt,msct

if (msct eq 'MON') then begin
  scaf1=1./(a(tasid('SIGM',b),ptnum)/a(tasid('MON',b),ptnum))^2.
  scaf2=1./(a(tasid('SIGM',b),ptnum+1)/a(tasid('MON',b),ptnum+1))^2.
  endif else begin
  scaf1=a(tasid(msct,b),ptnum)
  scaf2=a(tasid(msct,b),ptnum+1)
  endelse
toscaf=scaf1+scaf2
scaf1=scaf1/toscaf
scaf2=scaf2/toscaf               
d=fltarr(b.pars)
for i=0,b.pars-1 do begin
  d(i)=scaf1*a(i,ptnum)+scaf2*a(i,ptnum+1)
  if (i eq tasid('SIGI',b)) then $
    d(i)=sqrt((scaf1*a(i,ptnum))^2.+(scaf2*a(i,ptnum+1))^2.) else $
  if (i eq tasid('SIGM',b)) then $
    d(i)=sqrt((scaf1*a(i,ptnum))^2.+(scaf2*a(i,ptnum+1))^2.)
  if (a(i,ptnum) eq -9999.) then d(i)=a(i,ptnum+1)
  if (a(i,ptnum+1) eq -9999.) then d(i)=a(i,ptnum)
  endfor
c=fltarr(b.pars,npt-1)
if (ptnum gt 0) then c(*,0:ptnum-1)=a(*,0:ptnum-1)
c(*,ptnum)=d
if (ptnum lt npt-2) then c(*,ptnum+1:npt-2)=a(*,ptnum+2:npt-1)
return,c
end

function tassb,a,b,xind
;+
; NAME:         TASSB
; PURPOSE:
;               subtract two TASCOM data structures
; CATEGORY:
; CALLING SEQUENCE:
;               c=tassb(a,b,'x-index')
; INPUTS:
;               a,b- the two TASCOM data structures to be merged
;               'x-index'- the print parameter to use to index the
;                       data structures for rebining, ie if two points
;                       have the same value (within 5% of step size) of 
;                       x-index they will be subtracted
;                  e.g. c=tassb(a,b,'h') gives c=a-b in the intenisty column
; OPTIONAL INPUT PARAMETERS:
; KEYWORD PARAMETERS:
; OUPUTS:
;               a new data structure consisting of:
;                 d.points - integer - # of points in scan
;                 d.heads - integer - # of header parameters
;                 d.pars - integer - # of scan parameters
;                 d.text - string array - 2 lines of descriptive text
;                 d.hdef - string array - definition of header parameters
;                 d.hnum - float array - header parameters
;                 d.prp - string array - defintion of scan parameters
;                 d.data - 2D float array - scan data (parameter,point#)
; OPTIONAL OUTPUT PARAMETERS: 
; COMMON BLOCKS:
; SIDE EFFECTS:
;               if one of the TASCOM structures is missing a print parameter
;               field the other is used, only fields found in the first
;               structure (a) will be duplicated in the result (c).  Errors
;               are calculated RMS, only points found in both structures are
;               retained, fields other than i,sigi assume the values in a
; RESTRICTIONS:
; PROCEDURE:                                                               
; MODIFICATION HISTORY:
;               Thom Mason, Riso National Laboratory, 1993
;-
if (n_params() eq 2) then xind='H'
if (tasid(xind,a) eq -1) then xind=a.prp(0)
temp=a.text
if (strmid(a.text(0),0,2) eq strmid(b.text(0),0,2)) then $
  temp(0)=strcompress('('+a.text(0)+')'+'-('+strmid(b.text(0),3,100),/remove_all)+')' else $
  temp(0)=strcompress('('+a.text(0)+')'+'-'+b.text(0),/remove_all)
tol=abs(0.05*(a.data(tasid(xind,a),0)-a.data(tasid(xind,a),1)))
if (tasid('MON',a) eq -1) then msct='TIM' else $
  if (a.data(tasid('MON',a),0) le 0) then msct='TIM' else msct='MON'
msc=a.data(tasid(msct,a),0)
sortin=sort(a.data(tasid(xind,a),*))  ;first check to remove any duplicates in a
a.data(*,*)=a.data(*,sortin)
ptnum=0
npt=a.points
while (ptnum lt npt-1) do begin
  step=abs(a.data(tasid(xind,a),ptnum)-a.data(tasid(xind,a),ptnum+1))
  if (step le tol) then begin
    a.data=combine(a.data,ptnum,a,npt,msct)
    npt=npt-1
    endif else ptnum=ptnum+1
  endwhile
sortin=sort(b.data(tasid(xind,a),*))  ;now check to remove any duplicates in b
b.data(*,*)=b.data(*,sortin)
ptnum=0
npt=b.points
while (ptnum lt npt-1) do begin
  step=abs(b.data(tasid(xind,b),ptnum)-b.data(tasid(xind,b),ptnum+1))
  if (step le tol) then begin
    b.data=combine(b.data,ptnum,b,npt,msct)
    npt=npt-1
    endif else ptnum=ptnum+1
  endwhile
c=fltarr(a.pars,(a.points < b.points))     
ap=a.points
bp=b.points
numbo=0
for i=0,a.points-1 do begin
  bnum=where(abs(b.data(tasid(xind,b),*)-a.data(tasid(xind,a),i)) le tol)
  bnum=bnum(0)
  if (bnum ge 0) then begin
    c(*,numbo)=a.data(*,i)
    if (tasid('I',a) ge 0) then c(tasid('I',a),numbo)=(a.data(tasid('I',a),i)/a.data(tasid(msct,a),i)-$
      b.data(tasid('I',b),bnum)/b.data(tasid(msct,b),bnum))*msc
    if (tasid('SIGI',a) ge 0) then c(tasid('SIGI',a),numbo)=(((a.data(tasid('SIGI',a),i)/a.data(tasid(msct,a),i))^2.+$
      (b.data(tasid('SIGI',b),bnum)/b.data(tasid(msct,b),bnum))^2.)^0.5)*msc
    if (tasid('SIGM',a) ge 0) then c(tasid('SIGM',a),numbo)=msc*a.data(tasid('SIGM',a),i)/a.data(tasid(msct,a),i)
    if (tasid('MON',a) ge 0) then c(tasid('MON',a),numbo)=msc*a.data(tasid('MON',a),i)/a.data(tasid(msct,a),i)
    if (tasid('TIM',a) ge 0) then c(tasid('TIM',a),numbo)=msc*a.data(tasid('TIM',a),i)/a.data(tasid(msct,a),i)
    numbo=numbo+1
    endif
  end
new={points: numbo,heads: a.heads,pars: a.pars,text: temp,hdef: a.hdef,$
  hnum: a.hnum,prp: a.prp,data: c(*,0:numbo-1)}
return,new
end
