#!/software/perl5/bin/perl
#
  use Sys::Hostname;

  $host = hostname();

$arch_dir1 = "/disks/lv2131/apache/share/htdocs/MAG.Area/daily_data/raw/bnt/";

#$arch_dir1 = "/disks/de002/test_near_archive/";

 $total_files_proc = 0;


 if ( $host eq "spectrum" )
  {
   if($#ARGV < 2)
    {
      die " Usage: magview_end_to_end.pl <yyyy> <doy> <number of days>";
    }

   ($year,$doy,$num_days) = @ARGV;
#
   for ($count = 1; $count <= $num_days; $count++)
    {
     if ( ($doy == 366) && ($year == 1999) )
      {
       $year++;
       $doy = 1;
      }
     elsif ( ($doy == 367) && ($year == 2000) )
      {
       $year++;
       $doy =1;
      }

     $str_doy = sprintf("%03d",$doy);
     $line="$year\n$str_doy\n";
#
# Check to see if file is online. If not, issue warning
# and go on to the next day.
#
     if ($year <= 1999 )
      {
       $temp_year = $year - 1900;
      }
     else
      {
       $temp_year = $year - 2000;
      }

     $str_year = sprintf("%02d",$temp_year);
     $temp_in_file = $arch_dir1."b".$str_year.$str_doy.".dat.gz";

     if (-e $temp_in_file)
      {
#
#  Run magview_sdc
#
#>>> magview_sdc.pro should be run by the SDC
#
#       print (" ***** Running magview_sdc for $year $doy\n");
# 
#       open(IDL,"|/usr/local/bin/idl");
#
#       print IDL "\@magview_sdc.cmd\n";
#

#
#       print IDL "$line";
#       close(IDL);
#
#   Do the merging
#
       print (" ***** Running magview_merge for $year $doy\n");

       open(IDL,"|/usr/local/bin/idl");

       print IDL "\@magview_merge.cmd\n";

       print IDL "$line";
       close(IDL);

       $total_files_proc++;
      }   #if (-e $temp_file_in)
     else
      {
       print("**** $temp_in_file is *not* online. ***\n");
       print("**** Skipping to next day. ****\n");
      }
     $doy++;
   }  
 
   if ($total_files_proc > 0)
    {
#
#   Do the mrg_corrections
#
     print (" ***** Running mrg_corrections_v8 for all requested days.\n");

     system ("ls b9*.mrg b0*.mrg > mrg_files.lis");

     system ("mrg_corrections_v8");

# Now remove bhsk files: (BJA 2/13/00)

     system ("rm bhsk*.dat");
#
#   Do the xform
#
     print (" ***** Running  mag_eros_products(xform) for all requested days.\n");

     system ("ls b9*.fix b0*.fix > fix_files.lis");

     system ("mag_eros_products");

     print ("Gzipping output files.\n");

     system ("gzip *.mrg *.fix b*.dat");

     print (" **** Moving zipped output products to output subdirectory.\n");

     system ("mv *.gz output/");

     system ("rm mrg_files.lis fix_files.lis");
    }
   else
    {
     print("**** None of the requested days were online. ****\n");
    }
  }  # if ( $host eq "spectrum" )
 else
  {
    print "**** You're on $host. ****\n";
    print "**** You must be on SPECTRUM to run this script. ***\n";
  } 
function ICallheater_correction,b_tb,bx,by,bz,th0,th1,b_count,nh

b_tb = double(b_tb)
bx = float(bx)
by = float(by)
bz = float(bz)
th0 = double(th0)
th1 = double(th1)
b_count = long(b_count)
nh = long(nh)
ref = bytarr(8)

print,'in ICallheater_correction'
x = CALL_EXTERNAL('idlmerge_hp.sl','Callheater_correction',$
                   b_tb,bx,by,bz,th0,th1,b_count,nh,$
                   val=ref)
end
function ICallmerge_data,CTPMAG_time_diff,b_tb,b_count,ctp_ti,ctp_count,$
                         ri1,ri2,merged_elements_tb

CTPMAG_time_diff = double(CTPMAG_time_diff)
b_tb = double(b_tb)
b_count = long(b_count)
ctp_ti = double(ctp_ti)
ctp_count = long(ctp_count)
ri1 = long(ri1)
ri2 = long(ri2)
merged_elements_tb = long(merged_elements_tb)
ref = bytarr(8)

x = CALL_EXTERNAL('idlmerge_hp.sl','Callmerge_data',$
                   CTPMAG_time_diff,b_tb,b_count,$
                   ctp_ti,ctp_count,ri1,ri2,$
                   merged_elements_tb,val=ref)
end
FUNCTION a0fit,dt
    
    xfit = 0.
    a0tmp = 0.
    
    xfit=dt/250.
    if (xfit gt 10.) then begin
       xfit=10.
    endif
    if (xfit lt 1.6) then begin
       a0tmp=-0.25486/1.6*xfit
    endif else begin
       if (xfit lt 5.8231) then begin
          a0tmp=-0.37254 + 0.37333*xfit -0.24958*xfit^2 +0.043238*xfit^3 -0.0027206*xfit^4
       endif else begin
          if (xfit lt 7.8) then begin
             a0tmp=5.1021-2.0683*xfit+0.16777*xfit^2
          endif else begin
             a0tmp=0.71892-0.19576*xfit
          endelse
       endelse
    endelse
    RETURN, a0tmp
END
FUNCTION a1fit,dt
    
    xfit = 0.
    a0tmp = 0.
    
    xfit=dt/250.
    if (xfit gt 10) then begin
       xfit=10
    endif
    if (xfit lt 1.6) then begin
       a0tmp=0.091701/1.6*xfit
    endif else begin
       a0tmp=-0.023367 + 0.064802*xfit + 0.0075614*xfit^2 -0.002371*xfit^3 + 0.00026533*xfit^4
    endelse
    return, a0tmp
END
FUNCTION a2fit,dt
    
    xfit = 0.
    a0tmp = 0.
    
    xfit=dt/250.
    if (xfit gt 10) then begin
       xfit=10
    endif
    if (xfit lt 1.6) then begin
       a0tmp=0.26085/1.6*xfit
    endif else begin
       if (xfit lt 5.8231) then begin
          a0tmp=0.44343 -0.45007*xfit + 0.27928*xfit^2 -0.048141*xfit^3 + 0.0030137*xfit^4
       endif else begin
          if (xfit lt 7.8) then begin
             a0tmp= -5.0837 + 2.0626*xfit -0.16733*xfit^2
          endif else begin
             a0tmp=-0.69107 + 0.19221*xfit
          endelse
       endelse
    endelse
    return, a0tmp
END
FUNCTION a3fit,dt
    
    xfit = 0.
    a0tmp = 0.
    
    xfit=dt/250.
    if (xfit gt 10) then begin
       xfit=10
    endif
    if (xfit lt 1.8181) then begin
       a0tmp=0
    endif else begin
       if (xfit lt 6.8653) then begin
          a0tmp=0.16646 -0.10238*xfit +0.0079538*xfit^2 -0.0010071*xfit^3
       endif else begin
          a0tmp=-8.6052 +3.1004*xfit -0.36676*xfit^2 +0.012725*xfit^3
       endelse
    endelse
    return, a0tmp
END
;=============================================================================
;
;  Name: acquire_files.pro
;
;  Description: Acquire the input files requested.
;
;  Inputs: inDate - yyddd
;
;  Outputs: infiles - filenames
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro acquire_files,request_ident,inDate,infiles,no_file

  full_infiles = strarr(4)
  cp_commands = strarr(4)
  gunzip_commands = strarr(4)
  zip_ext = '.gz'
  ans = ' '

; Create input filenames.

  create_input_filenames,inDate,infiles

; Create full path filenames for first day

  create_fullpath_filenames,infiles,full_infiles

; Create copy commands.

  create_cp_commands,full_infiles,cp_commands

; Create gunzip commands.

  create_gunzip_commands,infiles,gunzip_commands

; Check if initial file exists

  find_result = findfile(full_infiles(0),count=file_cnt)

  if (file_cnt gt 0) then begin
   print,' Copying files.'
   copy_input_files,cp_commands
   print,' Gunzipping files.'
   gunzip_files,gunzip_commands
  endif else begin 
   no_file = 1
   if (request_ident lt 2) then begin
    print,' **** ',infiles(0)+zip_ext,' NOT on-line. ****'
    stop
   endif else begin
    print,'*** A gap has been detected, but the next day is not online. Continue? y/n'
;;;    read,ans
    ans = 'y'
    if ( (ans eq 'n') or (ans eq 'N') ) then begin
     stop
    endif else begin
     print,' **** OK GAP WILL *NOT* BE FILLED **** '
    endelse
   endelse
  endelse

end
;=============================================================================
;
;  Name: apply_heater_correction.pro
;
;  Description: Apply heater ON/OFF corrections to spacecraft magnetometer.
;
;  Inputs: bnt = magnetometer data
;          heater_log_times = heater ON/OFF times
;
;  Outputs: Structure containing the start and stop times.
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth      07/15/97   Original code
;  B. J. Anderson  10/01/97   Add X and Z corrections.
;
; $Id: apply_heater_correction.pro,v 1.3 1997/10/01 13:41:46 anderson Exp $
;
; $Log: apply_heater_correction.pro,v $
; Revision 1.3  1997/10/01 13:41:46  anderson
; Added X and Z corrections.
;
; Revision 1.2  1997/08/28 16:54:38  anderson
; Changed index to longint for large files.
;
; Revision 1.1  1997/07/17 17:20:15  tothba1
; Initial revision
;
;
;==========================================================================
;
    pro  apply_heater_correction,tb,bnt,heater_log_times
;
     i = 0L
     num_tb = 0L
     ton = 0.0D
     toff = 0.0D
     
     fix_flag = 1
     
     num_tb = n_elements(tb) - 1
;
     for i = 0L, num_tb do begin
      check_heater_on = -1
      check_heater_on = $
               where((heater_log_times.time1 lt tb(i)) and $
                     (heater_log_times.time2 gt tb(i)))
;
;     if (check_heater_on(0) gt -1) then begin
;      bnt(0,i) = bnt(0,i) + 0.1
;      bnt(1,i) = bnt(1,i) + 5.7
;      bnt(2,i) = bnt(2,i) - 0.77
;     endif

      if (check_heater_on(0) gt -1) then begin
         if (fix_flag eq 0) then begin
             bnt(0,i) = bnt(0,i) + 0.1
             bnt(1,i) = bnt(1,i) + 5.7
             bnt(2,i) = bnt(2,i) - 0.77
          endif else begin
             ton = heater_log_times.time1(check_heater_on)
             toff = heater_log_times.time2(check_heater_on)
         dby = by_htr(tb(i),ton,toff)
             bnt(0,i) = bnt(0,i) - 0.1/5.7*dby
             bnt(1,i) = bnt(1,i) - dby
             bnt(2,i) = bnt(2,i) + 0.77/5.7*dby
          endelse
      endif
      
     endfor
;
    end
;=============================================================================
;
;  Name: apply_mag_corrections.pro
;
;  Description: Apply corrections to spacecraft magnetometer.
;
;  Inputs: bnt = magnetometer data
;          bnt_raw = uncorrected magnetometer data
;          mag_corrections = corrections to apply
;          magc_strings = names of corrections
;          heater_log_times = heater ON/OFF times
;          mag_corrections_applied = whether or not corrections have
;                                    actually been applied
;
;  Outputs: Structure containing the start and stop times.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   07/15/97   Original code
;
; $Id: apply_mag_corrections.pro,v 1.1 1997/07/17 17:20:27 tothba1 Exp $
;
; $Log: apply_mag_corrections.pro,v $
; Revision 1.1  1997/07/17 17:20:27  tothba1
; Initial revision
;
;
;=============================================================================
;
    pro  apply_mag_corrections,tb,bnt,bnt_raw,mag_corrections,$
                                   magc_strings,heater_log_times,$
                                   mag_corrections_applied
;
    i = 0
        num_corrections = fix(total((mag_corrections gt -1)))
;
        if (num_corrections gt 0) then begin
     print,' ** Reverting to RAW bnt values first. **'
     bnt = bnt_raw
     for i = 0,(num_corrections - 1) do begin
       case i of
            0: begin
                 print,' ** Applying ',magc_strings[0],' **'
                 apply_heater_correction,tb,bnt,heater_log_times
               end
       endcase
     endfor
         mag_corrections_applied = 1;
        endif else begin
         print,' ** NO CORRECTIONS SELECTED -> NONE APPLIED **'
        endelse
;
    endFUNCTION by_htr,t,ton,toff

    xhtr = 0.
    dt_on = 0.
    dt_jmp = 0.
    btmp = 0.
    
    dt_on=toff-ton
    if (dt_on le 2650.0) then begin
       xhtr=2.*(t-ton)/dt_on - 1.
    endif else begin
       if ((t-ton) lt 2650.0) then begin
          xhtr=2.*(t-ton)/2650.0 - 1.0
       endif else begin
          xhtr=1.0
       endelse
    endelse
    btmp=a0fit(dt_on) + a1fit(dt_on)*xhtr + a2fit(dt_on)*xhtr^2 + a3fit(dt_on)*xhtr^3
    
    dt_jmp=dt_on
    if (dt_jmp gt 2500.0) then begin
       dt_jmp=2500.0
    endif
    btmp=btmp + 0.5*( -5.8882-3.7154e-05*dt_jmp -5.7037-3.0013e-05*dt_jmp )
    
    RETURN, btmp
END
;=============================================================================
;
;  Name: calc_seconds.pro
;
;  Description: Convert utc to seconds.
;
;  Inputs: day, hour, minute, seconds
;
;  Outputs: seconds
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro calc_seconds,day,hours,minutes,secs,tot_seconds
 
 tot_seconds = 0.0D

 tot_seconds = secs + (minutes * 60.0D) + (hours * 3600.0D) + (day * 86400.0D)

end
;=============================================================================
;
;  Name: calc_tb0.pro
;
;  Description: Calculate the the time starting point for the delta times.
;
;  Inputs: 
;
;  Outputs: seconds
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro calc_tb0,infiles,tb0

 temp_str = strarr(1)
 tday = 0
 thour = 0
 tmin = 0
 tsec = 0.
 btb = 0.
 bnt0 = 0.
 bnt1 = 0.
 bnt2 = 0.
 brng = 0 
 b_time = 0.0D
 ctp_time = 0.0D
 t0 = 0
 t1 = 0
 t2 = 0
 t3 = 0 
 t4 = 0 
 t5 = 0 
 t6 = 0 
 t7 = 0 
 t8 = 0 
 t9 = 0 

; Read in the first line of data from the b file
 openr, lun, infiles(0), /get_lun
 readf,lun,temp_str
 b_format = '(1x,i4,2(1x,i3),1x,f6.2,f15.2,3(f15.5),i2)'
 readf,lun,tday,thr,tmin,tsec,btb,bnt0,bnt1,bnt2,brng,format=b_format
 calc_seconds,tday,thr,tmin,tsec,b_time
 free_lun,lun

; Read in the first line of data from the ctp file
 openr, lun, infiles(1), /get_lun
 readf,lun,temp_str
 ctp_format = '(1x,i4,2(1x,i3),1x,f6.2,f15.2,10(1x,i6))'
 readf,lun,tday,thr,tmin,tsec,btb,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,format=ctp_format
 calc_seconds,tday,thr,tmin,tsec,ctp_time
 free_lun,lun

; Determine which time is earlier and use that as the starting point
 if (b_time le ctp_time) then begin
   tb0 = b_time
 endif else begin
   tb0 = ctp_time
 endelse

end
;=============================================================================
;
;  Name: copy_input_files.pro
;
;  Description: Issue system commands to copy input files to local directory.
;
;  Inputs: Copy commands
;
;  Outputs: None
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro copy_input_files,cp_commands

 for i = 0,3 do begin
  spawn,cp_commands(i)
 endfor

end
;=============================================================================
;
;  Name: create_cp_commands.pro
;
;  Description: Create copy commands to copy input data files to
;               the local directory.
;
;  Inputs: fill_infiles - fill path input filenames
;
;  Outputs: Copy commands
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro  create_cp_commands,full_infiles,cp_commands

 for i = 0,3 do begin
  cp_commands(i) = 'cp '+full_infiles(i)+' .'
 endfor

end
;=============================================================================
;
;  Name: create_fullpth_filenames.pro
;
;  Description: Create input data filenames with full directory path
;               and the gzip extension.
;
;  Inputs: infiles - input filenames
;
;  Outputs: Names of the input files with directory path and gzip extension.
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro  create_fullpath_filenames,infiles,full_infiles

 arch_dir1 = '/disks/lv2131/apache/share/htdocs/MAG.Area/daily_data/raw/bnt/';
 arch_dir2 = '/disks/lv2131/apache/share/htdocs/MAG.Area/daily_data/raw/ctp/';
 arch_dir3 = '/disks/lv2131/apache/share/htdocs/MAG.Area/daily_data/raw/bctp/'
 arch_dir4 = '/disks/lv2131/apache/share/htdocs/MAG.Area/daily_data/raw/bhsk/'
;test1
; arch_dir1 = '/homes/tothba1/sdc/test/';
; arch_dir2 = '/homes/tothba1/sdc/test/';
; arch_dir3 = '/homes/tothba1/sdc/test/';
;test2
; arch_dir1 = '/disks/de002/test_near_archive/';
; arch_dir2 = '/disks/de002/test_near_archive/';
; arch_dir3 = '/disks/de002/test_near_archive/';

 zip_ext = '.gz'

 full_infiles(0) = arch_dir1+infiles(0)+zip_ext
 full_infiles(1) = arch_dir2+infiles(1)+zip_ext
 full_infiles(2) = arch_dir3+infiles(2)+zip_ext
 full_infiles(3) = arch_dir4+infiles(3)+zip_ext

end
;=============================================================================
;
;  Name: create_gunzip_commands.pro
;
;  Description: Create gunzip commands to gunzip input data files 
;
;  Inputs: infiles - input filenames
;
;  Outputs: Copy commands
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro  create_gunzip_commands,infiles,gunzip_commands

 zip_ext = '.gz'

 for i = 0,3 do begin
  gunzip_commands(i) = 'gunzip '+infiles(i)+zip_ext
 endfor

end
;=============================================================================
;
;  Name: create_input_filenames.pro
;
;  Description: Create input data filenames.
;
;  Inputs: inDate - yyddd string
;
;  Outputs: Names of the input files.
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro  create_input_filenames,inDate,infiles

 infiles(0) = 'b'+inDate+'.dat'
 infiles(1) = 'ctp'+inDate+'.dat'
 infiles(2) = 'bctp'+inDate+'.dat'
 infiles(3) = 'bhsk'+inDate+'.dat'

end
;=============================================================================
;
;  Name: create_rm_commands.pro
;
;  Description: Create rm commands to delete input data files 
;
;  Inputs: infiles - input filenames
;
;  Outputs: Copy commands
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro create_rm_commands,infiles,rm_commands

 for i = 0,2 do begin
  rm_commands(i) = '\rm '+infiles(i)
 endfor

end
;=============================================================================
;
;  Name: fill_data_gap.pro
;
;  Description: Fill data gaps at end of a file.
;
;  Inputs: 
;
;  Outputs: 
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro fill_data_gap,ctp_tot_seconds,b_tot_seconds,infiles,tb0,b_day,b_hr,b_min,b_sec,b_tb,bnt,rng,b_count,ctp_day,ctp_hr,ctp_min,ctp_sec,ctp_ti,c,tube_temp,probe_temp,ctp_count,bctp_day,bctp_hr,bctp_min,bctp_sec,bctp_ti,b_ctp,bctp_count

 temp_str = strarr(1)
 tday = 0
 thour = 0
 tmin = 0
 tsec = 0.
 ttbi = 0.
 bnt0 = 0.
 bnt1 = 0.
 bnt2 = 0.
 brng = 0 
 t0 = 0
 t1 = 0
 t2 = 0
 t3 = 0 
 t4 = 0 
 t5 = 0 
 t6 = 0 
 t7 = 0 
 t8 = 0 
 t9 = 0 
 tube_t =0.
 probe_t = 0.
 bctp0 = 0.
 bctp1 = 0.
 bctp2 = 0.
 tti = 0.
 ttbi = 0.
 bctp_tot_seconds = 0.0D

 if (ctp_tot_seconds gt b_tot_seconds) then begin
  openr, lun, infiles(0), /get_lun
  readf,lun,temp_str
  b_format = '(1x,i4,2(1x,i3),1x,f6.2,f15.2,3(f15.5),i2)'
  while ( (NOT EOF(lun)) and (b_count le 599999) and $
          (ctp_tot_seconds gt b_tot_seconds) ) do begin
   readf,lun,tday,thr,tmin,tsec,ttbi,bnt0,bnt1,bnt2,brng,format=b_format
   calc_seconds,tday,thr,tmin,tsec,b_tot_seconds
   b_count = b_count + 1
   b_day(b_count) = tday
   b_hr(b_count) = thr
   b_min(b_count) = tmin
   b_sec(b_count) = tsec
   b_tb(b_count) = b_tot_seconds - tb0
   bnt(0,b_count) = bnt0
   bnt(1,b_count) = bnt1
   bnt(2,b_count) = bnt2
   rng(b_count) = brng
  endwhile

  if (b_count gt 599999) then begin
   print,' *** Busted an array boundary in ',infiles(0)
   print,' *** Counter = ',b_count
   stop
  endif

  free_lun,lun

 endif else begin

  openr, lun,infiles(1),/get_lun
  readf,lun,temp_str
  ctp_format = '(1x,i4,2(1x,i3),1x,f6.2,f15.2,10(1x,i6),2(1x,f6.2))'
  while ( (NOT EOF(lun)) and (ctp_count le 599999) and $
          (b_tot_seconds gt ctp_tot_seconds) ) do begin
   readf,lun,tday,thr,tmin,tsec,ttbi,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,tube_t,probe_t,format=ctp_format
   calc_seconds,tday,thr,tmin,tsec,ctp_tot_seconds
   ctp_count = ctp_count + 1
   ctp_day(ctp_count) = tday
   ctp_hr(ctp_count) = thr
   ctp_min(ctp_count) = tmin
   ctp_sec(ctp_count) = tsec
   ctp_ti(ctp_count) = ctp_tot_seconds - tb0
   c(0,ctp_count) = t0
   c(1,ctp_count) = t1
   c(2,ctp_count) = t2
   c(3,ctp_count) = t3
   c(4,ctp_count) = t4
   c(5,ctp_count) = t5
   c(6,ctp_count) = t6
   c(7,ctp_count) = t7
   c(8,ctp_count) = t8
   c(9,ctp_count) = t9
   tube_temp(ctp_count) = tube_t
   probe_temp(ctp_count) = probe_t
  endwhile

  if (ctp_count gt 599999) then begin
   print,' *** Busted an array boundary in ',infiles(1)
   print,' *** Counter = ',ctp_count
   stop
  endif

  free_lun,lun 

  openr, lun, infiles(2), /get_lun
  readf,lun,temp_str
  bctp_format = '(1x,i4,2(1x,i3),1x,f6.2,f15.2,3(f15.5))'
  while ( (NOT EOF(lun)) and (bctp_count le 599999) and $
          (bctp_count lt ctp_count) ) do begin
   readf,lun,tday,thr,tmin,tsec,tti,bctp0,bctp1,bctp2,format=bctp_format
   calc_seconds,tday,thr,tmin,tsec,bctp_tot_seconds
   bctp_count = bctp_count + 1
   bctp_day(bctp_count) = tday
   bctp_hr(bctp_count) = thr
   bctp_min(bctp_count) = tmin
   bctp_sec(bctp_count) = tsec
   bctp_ti(bctp_count) = bctp_tot_seconds - tb0
   b_ctp(0,bctp_count) = bctp0
   b_ctp(1,bctp_count) = bctp1
   b_ctp(2,bctp_count) = bctp2
  endwhile

  if (bctp_count gt 599999) then begin
   print,' *** Busted an array boundary in ',infiles(2)
   print,' *** Counter = ',bctp_count
   stop
  endif

  free_lun,lun

 endelse
end
;=============================================================================
;
;  Name: get_correction_data.pro
;
;  Description: Read in various data sets used to make corrections.
;
;  Inputs: tb0 = B field start time. This is subtracted from the
;                times read in get_heater_log.
;
;  Outputs: Structures containing the correction data.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   07/10/97   Original code
;
; $Id: get_correction_data.pro,v 1.3 1998/03/12 16:15:58 anderson Exp $
;
; $Log: get_correction_data.pro,v $
; Revision 1.3  1998/03/12 16:15:58  anderson
; Revised correction data per flyby results.
;
; Revision 1.2  1998/01/22 22:43:26  tothba1
; Minimized array
;
; Revision 1.1  1997/07/17 17:19:21  tothba1
; Initial revision
;
;
;=============================================================================
;
    pro get_correction_data,tb0,heater_log_times,nh
;
;
        heater_log_times = { heater_log, $
                time1 : dblarr(25000), $
                time2 : dblarr(25000) }
;
;--- Get the correction data for the heater.
;
    get_heater_log,tb0,heater_log_times,nh
;
    end
;=============================================================================
;
;  Name: get_heater_log.pro
;
;  Description: Read the file heater_log.dat. This file
;               was created by and is maintained by Brian
;               Anderson. Generally speaking it defines
;               when the heater was turned off and on.
;
;  Inputs: tb0 = B field start time. This is subtracted from the
;                heater times.
;
;  Outputs: Structure containing the start and stop times.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   07/10/97   Original code
;
; $Id: get_heater_log.pro,v 1.2 1998/03/12 16:16:42 anderson Exp $
;
; $Log: get_heater_log.pro,v $
; Revision 1.2  1998/03/12 16:16:42  anderson
; Revised per flyby results.
;
; Revision 1.1  1998/01/22 22:42:51  tothba1
; Initial revision
;
; Revision 1.1  1997/07/17 17:19:40  tothba1
; Initial revision
;
;
;
;=============================================================================
;
    pro get_heater_log,tb0,heater_log_times,nh
;
;
    i = 0L
    temp_time1 = 0.D0
    temp_time2 = 0.D0
    heater_log_times = { heater_log, $
                time1 : dblarr(25000), $
                time2 : dblarr(25000) }
;
;--- Get a logical unit number and open heater_log.dat
;
    get_lun,lun
    openr,lun,'heater_log.dat'
;
;--- Read in the heater times.
;
    while (not eof(lun)) do begin
;
     readf,lun,format='(2(2x,f13.2))',temp_time1, $
                                          temp_time2
         temp_time1 = temp_time1 - tb0
         temp_time2 = temp_time2 - tb0
;
         if( (temp_time1 ge 0.0) or (temp_time2 ge 0.0) )then begin
           heater_log_times.time1(i) = temp_time1
           heater_log_times.time2(i) = temp_time2
       i = i + 1
         endif
;
    endwhile
;
        nh = i - 1
    close,lun
;
    end


;=============================================================================
;
;  Name: gunzip_files.pro
;
;  Description: Issue system commands to gunzip input files.
;
;  Inputs: Gunzip commands
;
;  Outputs: None
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro gunzip_files,gunzip_commands

 for i = 0,3 do begin
  spawn,gunzip_commands(i)
 endfor

end
;*****************************************************************************
;
;  Name: magview_merge.pro
;
;  Description: Merge the mag and ctp data.
;
;  Inputs: year and doy
;
;  Outputs: Merged data files.
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Brian Anderson  1996-97    Original code
;  Bruce Toth       4/1999    Reworked to handle data gaps
;
;-----------------------------------------------------------------------;

@ICallmerge_data.pro
@ICallheater_correction.pro

  infiles = strarr(4)

  heater_log_times = { heater_log, $
                  time1 : dblarr(25000), $
                  time2 : dblarr(25000) }

  nh = 0L    ;number of heater log entries
  request_ident = 1
  no_file = 0
  b_header = strarr(1)
  b_day = intarr(600000)
  b_hr = intarr(600000)
  b_min = intarr(600000)
  b_sec = fltarr(600000)
  tb0 = 0.0D
  b_tb = dblarr(600000)
  bnt = fltarr(3,600000)
  rng = intarr(600000)
  b_count = 0L
  ctp_header = string(1)
  ctp_day = intarr(600000)
  ctp_hr = intarr(600000)
  ctp_min = intarr(600000)
  ctp_sec = fltarr(600000)
  ctp_ti = dblarr(600000)
  c = intarr(10,600000)
  tube_temp = fltarr(600000)
  probe_temp = fltarr(600000)
  ctp_count = 0L
  bctp_header = strarr(1)
  bctp_day = intarr(600000)
  bctp_hr = intarr(600000)
  bctp_min = intarr(600000)
  bctp_sec = fltarr(600000)
  bctp_ti = dblarr(600000)
  b_ctp = fltarr(3,600000)
  bctp_count = 0L
  rm_commands = strarr(3)
  b_tot_seconds = 0.0D
  ctp_tot_seconds = 0.0D

  outstr=['To screen','Post-script File']

  inDate = ' '
  in2YR = ' '
  inYR = ' '
  inDOY = ' '
;  print,'Enter date of data to read: YYYY,DDD'
  read,inYR,inDOY
  
  intDOY1 = FIX(inDOY)
  intDOY2 = intDOY1 + 1
  inDOY2 = strcompress(string(intDOY2),/REMOVE_ALL)
  intYR=FIX(inYR)
  sto_intYR = intYR
;
; Construct date string for intial file.
;
  if (intYR lt 2000) then begin
     intYR=intYR-1900
     inDate = strcompress(STRING(format='(i,a3)',intYR,inDOY),/REMOVE_ALL)
  endif else begin  ; date later than 2000
     intYR=intYR-2000
     inDate = STRING(format='(i2.2,a3)',intYR,inDOY)
  endelse
;
; Construct an additional date string in case more records are
; required from a second file.
;
  if ( intDOY2 le 365) then begin
    inDate2 = STRING(format='(i2.2,i3.3)',intYR,intDOY2)  ;inDOY2)
  endif else begin  
   if ( (sto_intYR + 1) lt 2000) then begin  
    intYR2 = (sto_intYR + 1) - 1900
   endif else begin
    intYR2 = (sto_intYR + 1) - 2000
   endelse
   inDOY2 = '001'
   inDate2 = STRING(format='(i2.2,a3)',intYR2,inDOY2)
  endelse

; Copy and gunzip requested input files in local directory
  acquire_files,request_ident,inDate,infiles,no_file

; Create the rm commands for later use.
  create_rm_commands,infiles,rm_commands

; Determine tb0
  calc_tb0,infiles,tb0

; Read in the byyddd.dat data
  read_b_data,infiles(0),tb0,b_header,b_day,b_hr,b_min,b_sec,b_tb,bnt,rng,b_count

; Read in the ctpyyddd.dat data
  read_ctp_data,infiles(1),tb0,ctp_header,ctp_day,ctp_hr,ctp_min,ctp_sec,ctp_ti,c,tube_temp,probe_temp,ctp_count

; Read in the bctpyyddd.dat data.
  read_bctp_data,infiles(2),tb0,bctp_header,bctp_day,bctp_hr,bctp_min,bctp_sec,bctp_ti,b_ctp,bctp_count

; We can delete the input data files now.
  rm_files,rm_commands

  calc_seconds,b_day(b_count),b_hr(b_count),b_min(b_count),b_sec(b_count),$
               b_tot_seconds
  calc_seconds,ctp_day(ctp_count),ctp_hr(ctp_count),ctp_min(ctp_count),$
               ctp_sec(ctp_count),ctp_tot_seconds


; Fill in any end of file data gaps.
  if (ctp_tot_seconds gt b_tot_seconds) then begin
   print,'Attempting to fill data gap....ctp time > b time'
   request_ident = 2
   acquire_files,request_ident,inDate2,infiles,no_file
   if ( no_file lt 1 ) then begin
    create_rm_commands,infiles,rm_commands
    print,'Attempting to fill data gap.'
    fill_data_gap,ctp_tot_seconds,b_tot_seconds,infiles,tb0,b_day,b_hr,b_min,b_sec,b_tb,bnt,rng,b_count,ctp_day,ctp_hr,ctp_min,ctp_sec,ctp_ti,c,tube_temp,probe_temp,ctp_count,bctp_day,bctp_hr,bctp_min,bctp_sec,bctp_ti,b_ctp,bctp_count
    rm_files,rm_commands
   endif
  endif else begin
   if (b_tot_seconds gt ctp_tot_seconds) then begin
    print,'Attempting to fill data gap....b time > ctp time'
    request_ident = 2
    acquire_files,request_ident,inDate2,infiles,no_file
    if ( no_file lt 1 ) then begin
     create_rm_commands,infiles,rm_commands
     fill_data_gap,ctp_tot_seconds,b_tot_seconds,infiles,tb0,b_day,b_hr,b_min,b_sec,b_tb,bnt,rng,b_count,ctp_day,ctp_hr,ctp_min,ctp_sec,ctp_ti,c,tube_temp,probe_temp,ctp_count,bctp_day,bctp_hr,bctp_min,bctp_sec,bctp_ti,b_ctp,bctp_count
     rm_files,rm_commands
    endif
   endif
  endelse

; Resize input arrays
  b_tb = b_tb(0:b_count)
  bnt = bnt(*,0:b_count)
  ctp_ti = ctp_ti(0:ctp_count)
  bctp_ti = bctp_ti(0:bctp_count)

  hr0 = -2 ;0
  hr1 = 26 ;24

  tstart = double(3600.0*hr0)
  tstop = double(3600.0*hr1)
;
;--- Get various kinds of correction data.
;
  get_correction_data,tb0,heater_log_times,nh
 
  mag_corrections = intarr(1)  ; mag corrections
  mag_corrections(*) = -1      ; none selected
  mag_corrections(0) = 0       ; default-> heater correction ON
  mag_corrections_applied  = 1 ; because we start with heater correction ON<
  magc_strings=['Heater Jump Corrections.']
;
;--- We apply heater corrections by default because Brian
;    said so.

  bnt_raw = bnt
;  apply_mag_corrections,b_tb,bnt,bnt_raw,mag_corrections,$
;                        magc_strings,heater_log_times

   th0 = heater_log_times.time1(0:nh)
   th1 = heater_log_times.time2(0:nh)
   bx = bnt(0,*)
   by = bnt(1,*)
   bz = bnt(2,*)
;
   y = ICallheater_correction(b_tb,bx,by,bz,th0,th1,b_count,nh)

   bnt(0,*) = bx
   bnt(1,*) = by
   bnt(2,*) = bz
;
;--- Let's merge the mag and the CTP data - this must be done after the
;    heater corrections are applied.
;
  CTPMAG_time_diff = ' '
;;  print,'Enter Maximum CTP-MAG time difference in seconds:'
;;  read,CTPMAG_time_diff
;;  flt_CTPMAG_time_diff = float(CTPMAG_time_diff)

  CTPMAG_time_diff = 50.0D
;  coin_index,CTPMAG_time_diff,b_tb,ctp_ti,ri1,ri2

   ri1  = intarr(600000)
   ri2  = intarr(600000)
   merged_elements_tb = 0L

   x = ICallmerge_data(CTPMAG_time_diff,b_tb,b_count,ctp_ti,ctp_count,$
                       ri1,ri2,merged_elements_tb)
;
;--- Now we create the merged arrays
;
;  merged_elements_tb = n_elements(ri1)
;  merged_elements_ti = n_elements(ri2)

  ri1 = ri1(0:merged_elements_tb-1)
  ri2 = ri2(0:merged_elements_tb-1)
  merged_elements_ti = merged_elements_tb
  bnt_merged = fltarr(3,merged_elements_tb)
  c_merged = intarr(10, merged_elements_ti)
  b_ctp_merged = fltarr(3,merged_elements_ti)

  tb_m = b_tb(ri1)
  ti_m = ctp_ti(ri2)

  b_day = b_day(ri1)
  b_hr = b_hr(ri1)
  b_min = b_min(ri1)
  b_sec = b_sec(ri1)
  rng_merged = rng(ri1)
  ctp_day = ctp_day(ri2)
  ctp_hr = ctp_hr(ri2)
  ctp_min = ctp_min(ri2)
  ctp_sec = ctp_sec(ri2)
  bctp_day = bctp_day(ri2)
  bctp_hr = bctp_hr(ri2)
  bctp_min = bctp_min(ri2)
  bctp_sec = bctp_sec(ri2)
  bctp_ti = bctp_ti(ri2)

  bavg=fltarr(3)

  for k = 0L,merged_elements_tb-1 do begin
     navg=1
     bavg(*)=bnt(0:2,ri1(k))
     
     ia=1
avgloop1:
     if ((ia+ri1(k)) lt b_count ) then begin
        if ( b_tb(ri1(k)+ia) le (b_tb(ri1(k))+0.5) ) then begin
       bavg(*)=bavg(*)+bnt(0:2,ri1(k)+ia)
       navg=navg+1
       ia=ia+1
       goto,avgloop1
    endif
     endif
     
     ia=-1
avgloop2:
     if ((ia+ri1(k)) gt 0 ) then begin
        if ( b_tb(ri1(k)+ia) ge (b_tb(ri1(k))-0.5) ) then begin
       bavg(*)=bavg(*)+bnt(0:2,ri1(k)+ia)
       navg=navg+1
       ia=ia-1
       goto,avgloop2
    endif
     endif
     
     bavg=bavg/float(navg)
     bnt_merged(*,k)=bavg(0:2)
     
  endfor

  for ii = 0,9 do begin
   c_merged(ii,*) = c(ii,ri2)
  endfor

  for ii = 0,2 do begin
   b_ctp_merged(ii,*) = b_ctp(ii,ri2)
  endfor

;
; Write output files.   
;
  ASCII_file='b'+inDate+'.mrg'
  openw, 1, ASCII_file
  printf,1,format='("Day Hr Mn Sec   t-t0 (sec)   BX, BY, BZ, Range")'
  i=long(0)
  while (i le (merged_elements_tb - 1)) do begin
    printf,1,format='(1x,i4,2(1x,i3),1x,f6.2,f15.2,3(f15.5),i2)',$
              b_day(i),b_hr(i),b_min(i),b_sec(i),tb_m(i),$
              bnt_merged(0,i),bnt_merged(1,i),bnt_merged(2,i),rng_merged(i)
              i=i+long(1)
  endwhile
  close,1

  ASCII_file='ctp'+inDate+'.mrg'
  openw, 1, ASCII_file
  printf,1,ctp_header
        
  i=long(0)
  while (i le (merged_elements_ti - 1)) do begin
    printf,1,format='(1x,i4,2(1x,i3),1x,f6.2,f15.2,$)',$
    ctp_day(i),ctp_hr(i),ctp_min(i),ctp_sec(i),ti_m(i)
    j=0
    while (j lt 10) do begin
      if (j le 9) then printf,1,format='($,1x,i6)',c_merged(j,i)
;      if (j lt 9) then printf,1,format='($,1x,i6)',c_merged(j,i)
;;;      if (j eq 9) then printf,1,format='(1x,i6)',c_merged(j,i)
      j=j+1
    endwhile
    printf,1,format='($,1x,f7.2)',tube_temp(i)
    printf,1,format='(1x,f7.2)',probe_temp(i)
    i=i+long(1)
  endwhile
  close,1

  ASCII_file='bctp'+inDate+'.mrg'
  openw, 1, ASCII_file
  printf,1,format='("Day Hr Mn Sec   t-t0 (sec),BX_ctp,BY_ctp,BZ_ctp")'
  i=long(0)
  while (i le (merged_elements_ti - 1)) do begin
    printf,1,format='(1x,i4,2(1x,i3),1x,f6.2,f15.2,3(f15.5))',$
             bctp_day(i),bctp_hr(i),bctp_min(i),bctp_sec(i),$
             bctp_ti(i),b_ctp_merged(0,i),b_ctp_merged(1,i),b_ctp_merged(2,i)
    i=i+long(1)
  endwhile
  close,1

end  ;Main
;=============================================================================
;
;  Name: read_b_data.pro
;
;  Description: Read the byyddd.dat file.
;
;  Inputs: filename
;
;  Outputs: Data arrays
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro read_b_data,filename,tb0,b_header,b_day,b_hr,b_min,b_sec,b_tb,bnt,rng,b_count

 i = 0L
 bday = 0
 bhour = 0
 bmin = 0
 bsec = 0.
 btb = 0.
 calc_time = 0.0D
 bnt0 = 0.
 bnt1 = 0.
 bnt2 = 0.
 brng = 0 

 openr, lun, filename, /get_lun

 readf,lun,b_header
 b_format = '(1x,i4,2(1x,i3),1x,f6.2,f15.2,3(f15.5),i2)'
 while ( (NOT EOF(lun)) and (i le 599999) ) do begin
  readf,lun,bday,bhr,bmin,bsec,btb,bnt0,bnt1,bnt2,brng,format=b_format
  calc_seconds,bday,bhr,bmin,bsec,calc_time
  b_day(i) = bday
  b_hr(i) = bhr
  b_min(i) = bmin
  b_sec(i) = bsec
  b_tb(i) = calc_time - tb0
  bnt(0,i) = bnt0
  bnt(1,i) = bnt1
  bnt(2,i) = bnt2
  rng(i) = brng
  i = i + 1
 endwhile

 if (i gt 599999) then begin
  print,' *** Busted an array boundary in ',filename
  print,' *** Counter = ',i
  stop
 endif

 b_count = i - 1 

 free_lun,lun

end
;=============================================================================
;
;  Name: read_bctp_data.pro
;
;  Description: Read the bctpyyddd.dat file.
;
;  Inputs: filename
;
;  Outputs: Data arrays
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro read_bctp_data,filename,tb0,bctp_header,bctp_day,bctp_hr,bctp_min,bctp_sec,bctp_ti,b_ctp,bctp_count

 i = 0L
 calc_time = 0.0D
 bctpday = 0
 bctphour = 0
 bctpmin = 0
 bctpsec = 0.
 bctpti = 0.
 bctp0 = 0.
 bctp1 = 0.
 bctp2 = 0.
 
 openr, lun, filename, /get_lun

 readf,lun,bctp_header
 bctp_format = '(1x,i4,2(1x,i3),1x,f6.2,f15.2,3(f15.5))'
 while ( (NOT EOF(lun)) and (i le 599999) ) do begin
  readf,lun,bctpday,bctphr,bctpmin,bctpsec,bctpti,bctp0,bctp1,bctp2,format=bctp_format
  calc_seconds,bctpday,bctphr,bctpmin,bctpsec,calc_time
  bctp_day(i) = bctpday
  bctp_hr(i) = bctphr
  bctp_min(i) = bctpmin
  bctp_sec(i) = bctpsec
  bctp_ti(i) = calc_time - tb0
  b_ctp(0,i) = bctp0
  b_ctp(1,i) = bctp1
  b_ctp(2,i) = bctp2
  i = i + 1
 endwhile

 if (i gt 599999) then begin
  print,' *** Busted an array boundary in ',filename
  print,' *** Counter = ',i
  stop
 endif

 bctp_count = i - 1

 free_lun,lun

end
;=============================================================================
;
;  Name: read_ctp_data.pro
;
;  Description: Read the ctpyyddd.dat file.
;
;  Inputs: filename
;
;  Outputs: Data arrays
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro read_ctp_data,filename,tb0,ctp_header,ctp_day,ctp_hr,ctp_min,ctp_sec,ctp_ti,c,tube_temp,probe_temp,ctp_count

 i = 0L
 ctpday = 0
 ctphour = 0
 ctpmin = 0
 ctpsec = 0.
 ctpti = 0.
 t0 = 0
 t1 = 0
 t2 = 0
 t3 = 0 
 t4 = 0 
 t5 = 0 
 t6 = 0 
 t7 = 0 
 t8 = 0 
 t9 = 0 
 tube_t = 0.
 probe_t = 0.
 calc_time = 0.0D

 openr, lun, filename, /get_lun

 readf,lun,ctp_header
 ctp_format = '(1x,i4,2(1x,i3),1x,f6.2,f15.2,10(1x,i6),2(1x,f7.2))'
 while ( (NOT EOF(lun)) and (i le 599999) ) do begin
  readf,lun,ctpday,ctphr,ctpmin,ctpsec,ctpti,t0,t1,t2,t3,t4,t5,t6,$
            t7,t8,t9,tube_t,probe_t,format=ctp_format
  calc_seconds,ctpday,ctphr,ctpmin,ctpsec,calc_time
  ctp_day(i) = ctpday
  ctp_hr(i) = ctphr
  ctp_min(i) = ctpmin
  ctp_sec(i) = ctpsec
  ctp_ti(i) = calc_time - tb0
  c(0,i) = t0
  c(1,i) = t1
  c(2,i) = t2
  c(3,i) = t3
  c(4,i) = t4
  c(5,i) = t5
  c(6,i) = t6
  c(7,i) = t7
  c(8,i) = t8
  c(9,i) = t9
  tube_temp(i) = tube_t
  probe_temp(i) = probe_t
  i = i + 1
 endwhile

 if (i gt 599999) then begin
  print,' *** Busted an array boundary in ',filename
  print,' *** Counter = ',i
  stop
 endif

 ctp_count = i - 1 

 free_lun,lun

end
;=============================================================================
;
;  Name: rm_files.pro
;
;  Description: Issue system commands to delete input files.
;
;  Inputs: rm commands
;
;  Outputs: None
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/99     Original code
;
; $Id:$
;
; $Log: $
;
;==========================================================================
;
pro rm_files,rm_commands

 for i = 0,2 do begin
  spawn,rm_commands(i)
 endfor

end
#include <stdio.h>

void heater_correction(double tb[],float bx[],float by[],float bz[],
                       double th0[], double th1[], int nb, int nh);

long Callheater_correction(int argc, void *argv[])
{
 double *tb;
 float *bx;
 float *by;
 float *bz;
 double *th0;
 double *th1;
 int nb;
 int nh;

 tb = (double *)argv[0];
 bx = (float *)argv[1];
 by = (float *)argv[2];
 bz = (float *)argv[3];
 th0 = (double *)argv[4];
 th1 = (double *)argv[5];
 nb = *(int *)argv[6];
 nh = *(int *)argv[7];

 heater_correction(tb,bx,by,bz,th0,th1,nb,nh);

 return((long)1);
}    
#include <stdio.h>

void merge_data(double dt_max,double tb[],long int na,
                double tc[],long int nb,int ri1[],
                int ri2[],long int *k);

long Callmerge_data(int argc, void *argv[])
{
 double dt_max; 
 double *tb;
 long int na;
 double *tc;
 long int nb;
 int *ri1;
 int *ri2;
 long int *k;

 dt_max = *(double *)argv[0];
 tb = (double *)argv[1];
 na = *(long int *)argv[2];
 tc = (double *)argv[3];
 nb = *(long int *)argv[4];
 ri1 = (int *)argv[5];
 ri2 = (int *)argv[6];
 k = (long int*)argv[7];

 merge_data(dt_max,tb,na,tc,nb,ri1,ri2,k);

 return((long)1);
}    
#include <math.h>

float a0fit(double dt)
{
 double xfit;
 float a0tmp;

 xfit = 0.;
 a0tmp = 0.;

 xfit = dt/250.;

 if (xfit > 10.) xfit = 10.;

 if (xfit < 1.6) 
  a0tmp = -0.25486/1.6*xfit;
 else
  {
   if (xfit < 5.8231) 
    a0tmp = -0.37254 + 0.37333*xfit - 0.24958*pow(xfit,2) + 0.043238*pow(xfit,3) - 0.0027206*pow(xfit,4);
   else
    {
     if (xfit < 7.8) 
      a0tmp = 5.1021 - 2.0683*xfit + 0.16777*pow(xfit,2);
     else
      a0tmp = 0.71892 - 0.19576*xfit;
    }
  }
 return a0tmp;
}
#include <math.h>

float a1fit(double dt)
{
 double xfit;
 float a0tmp;

 xfit = 0.;
 a0tmp = 0.;

 xfit = dt/250.;

 if (xfit >10) xfit = 10;

 if (xfit < 1.6)
  a0tmp = 0.091701/1.6*xfit;
 else
  a0tmp = -0.023367 + 0.064802*xfit + 0.0075614*pow(xfit,2) - 0.002371*pow(xfit,3) + 0.00026533*pow(xfit,4);
 
 return a0tmp;
}
#include <math.h>

float a2fit(double dt)
{
  double xfit;
  float a0tmp;

  xfit = 0.;
  a0tmp = 0.;
    
  xfit = dt/250.;

  if (xfit > 10) xfit = 10;

  if (xfit < 1.6) 
   a0tmp = 0.26085/1.6*xfit;
  else
   {
    if (xfit < 5.8231)
     a0tmp=0.44343 - 0.45007*xfit + 0.27928*pow(xfit,2) - 0.048141*pow(xfit,3) + 0.0030137*pow(xfit,4);
    else
     {
      if (xfit < 7.8)
       a0tmp = -5.0837 + 2.0626*xfit - 0.16733*pow(xfit,2);
      else
       a0tmp = -0.69107 + 0.19221*xfit;
     }
   }
  return a0tmp;
}
#include <math.h>

float a3fit(double dt)
{
 double xfit;
 float a0tmp;

 xfit = 0.;
 a0tmp = 0.;

 xfit = dt / 250.;

 if (xfit > 10) xfit=10;

 if (xfit < 1.8181)
   a0tmp=0;
 else
  {
   if (xfit < 6.8653)
    a0tmp= 0.16646 - 0.10238*xfit + 0.0079538*pow(xfit,2) -0.0010071*pow(xfit,3);
   else
    a0tmp = -8.6052 + 3.1004*xfit - 0.36676*pow(xfit,2) + 0.012725*pow(xfit,3);
  }
 return a0tmp;
}
/* ===========================================================================
;
;  Name: build_bctp_ctp_file_names
;
;  Description: Build the bctp & ctp file names based on the byyddd.fix files.
;
;  Inputs: List of byyddd.fix files.
;
;  Outputs: List of bctp* files.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/15/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>
#include <string.h>

void build_bctp_ctp_file_names(int num_b_files,char b_input_files[100][50],
                               char bctp_input_files[100][50],
                               char ctp_input_files[100][50])
{
 int i = 0;
 char bctp_root[5] = "bctp";
 char ctp_root[4] = "ctp";
 char bctp_ctp_ext[5] = ".mrg";

 for (i = 0; i <= num_b_files; ++i)
  {
/* Construct the bctpyyddd.mrg input file. */

   strcpy(bctp_input_files[i],"             ");
   strcpy(bctp_input_files[i],bctp_root);
   bctp_input_files[i][4] = b_input_files[i][1];
   bctp_input_files[i][5] = b_input_files[i][2];
   bctp_input_files[i][6] = b_input_files[i][3];
   bctp_input_files[i][7] = b_input_files[i][4];
   bctp_input_files[i][8] = b_input_files[i][5];

   bctp_input_files[i][9] = bctp_ctp_ext[0];
   bctp_input_files[i][10] = bctp_ctp_ext[1];
   bctp_input_files[i][11] = bctp_ctp_ext[2];
   bctp_input_files[i][12] = bctp_ctp_ext[3];

/* Construct the ctpyyddd.mrg input file. */

   strcpy(ctp_input_files[i],"            ");
   strcpy(ctp_input_files[i],ctp_root);
   ctp_input_files[i][3] = b_input_files[i][1];
   ctp_input_files[i][4] = b_input_files[i][2];
   ctp_input_files[i][5] = b_input_files[i][3];
   ctp_input_files[i][6] = b_input_files[i][4];
   ctp_input_files[i][7] = b_input_files[i][5];

   ctp_input_files[i][8] = bctp_ctp_ext[0];
   ctp_input_files[i][9] = bctp_ctp_ext[1];
   ctp_input_files[i][10] = bctp_ctp_ext[2];
   ctp_input_files[i][11] = bctp_ctp_ext[3];

  }
}

/* ===========================================================================
;
;  Name: build_output_file_name
;
;  Description: Construct the name of output files.
;
;  Inputs: year, doy
;
;  Outputs: Name of the output files.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/15/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <string.h>

void build_output_file_names(int year, int doy, char output_file_names[4][50])
{
 char intro_string[3] = "b_";
 char file_ident[4][8] = {"SC_","J2000_","EBF_","NSO_"};
 char temp_year[5];
 char temp_doy[4];
 char file_ext[5] = ".dat";
 int i;

 if ( year <= 1999 )
  year = year - 1900;
 else
  year = year - 2000;

 sprintf(temp_year,"%02d",year);
 sprintf(temp_doy,"%03d",doy);

 for (i = 0; i <= 3; i++)
  {
   strcpy(output_file_names[i],intro_string);
   strcat(output_file_names[i],file_ident[i]);
   strcat(output_file_names[i],temp_year);
   strcat(output_file_names[i],temp_doy);
   strcat(output_file_names[i],file_ext);
  }
}
#include <math.h>

float a0fit(double dt);
float a1fit(double dt);
float a2fit(double dt);
float a3fit(double dt);

float by_htr(double t, double ton, double toff)
{
 double dt_on;
 double xhtr,dt_jmp,btmp;

 xhtr = 0.;
 dt_on = 0.;
 dt_jmp = 0.;
 btmp = 0.;

 dt_on = toff - ton;

 if (dt_on <= 2650.0)
  xhtr = 2.*(t-ton)/dt_on - 1.;
 else
  {
   if ((t-ton) < 2650.0)
    xhtr = 2.*(t-ton)/2650.0 - 1.0;
   else
    xhtr = 1.0;
  }

 btmp = a0fit(dt_on) + a1fit(dt_on)*xhtr + a2fit(dt_on)*pow(xhtr,2) + a3fit(dt_on)*pow(xhtr,3);
    
 dt_jmp = dt_on;

 if (dt_jmp > 2500.0)
  dt_jmp=2500.0;
    
 btmp = btmp + 0.5*(-5.8882 - 3.7154*pow(10,-5)*dt_jmp - 5.7037 - 3.0013*pow(10,-5)*dt_jmp);

 return btmp;
}
/* ===========================================================================
;
;  Name: calc_MET_time
;
;  Description: Calculate the MET time
;
;  Inputs: doy, hour, minutes, seconds
;
;  Outputs: MET
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   01/13/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */



void calc_MET_time(int dday, int hhr, int mmin, double ssec,double *cbMET)
{
 *cbMET = (dday * 86400.0) + (hhr * 3600.0) + (mmin * 60.0) + ssec;
}
/* ===========================================================================
;
;  Name: calc_load_times
;
;  Description: Calculate and load UTC and MET times.
;
;  Inputs: array counter, b_mrg and bctp_mrg data structures
;
;  Outputs: times data structure
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   12/17/97   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>

#include "spice_id.h"
#include "spice2.h"

#include "mag_eros_products_data_types.h"
#include "mag_eros_products_prototypes.h"

/* ---------------------------------------------------------------------- */
void calc_load_times(long int array_count, struct b_data b_mrg,
                     struct bctp_data bctp_mrg, struct time_data *proc_time,
                     struct time_errors *proc_time_errors)
{
 long int i;
 double b_MET, bctp_MET;
 int temp_day,temp_hr,temp_min;
 double temp_sec;
 double utc_time[3]; /* returns year, day, and seconds of day */
 int utc_hr, utc_min;
 float utc_sec;
 int errst;
 char utcstr[120];
 char *cptr;
 float fp_temp;
 int met_err;

 for (i = 0; i <= array_count; ++i)
  {
   /* Calc MET times and the delta MET times. */

   temp_day = b_mrg.day[i];
   temp_hr = b_mrg.hr[i];
   temp_min = b_mrg.min[i];
   temp_sec = b_mrg.sec[i]; 
   calc_MET_time(temp_day,temp_hr,temp_min,temp_sec,&b_MET);

   temp_day = bctp_mrg.day[i];
   temp_hr = bctp_mrg.hr[i];
   temp_min = bctp_mrg.min[i];
   temp_sec = bctp_mrg.sec[i]; 
   calc_MET_time(temp_day,temp_hr,temp_min,temp_sec,&bctp_MET); 

   proc_time->MET[i] = b_MET;
   proc_time->MET_diff[i] = fabs(b_MET - bctp_MET);

   /* Calc UTC time */

   met_err = metutc(NEAR, b_MET,utcstr,&errst);

   proc_time_errors->metutc_error[i] = met_err;
   cptr = strtok(utcstr,"-");
   proc_time->utc_year[i] = atoi(cptr);
   cptr = strtok(NULL," ");
   proc_time->utc_doy[i] = atoi(cptr);
   cptr = strtok(NULL," ");
   cptr = strtok(NULL,":");
   proc_time->utc_hr[i] = atoi(cptr);
   cptr = strtok(NULL,":");
   proc_time->utc_min[i] = atoi(cptr);
   cptr = strtok(NULL,":");
   proc_time->utc_sec[i] = (float) atof(cptr);
  }
}
/* ===========================================================================
;
;  Name: check_array_counters
;
;  Description: Make sure that 3 arrays have the same number of elements.
;               If there is a discrepancy, flag it and exit.
;
;  Inputs: Three array counters and file/array names.
;
;  Outputs: None.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/15/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>

void check_array_counters(long int array_counters[3], char filename_1[50], 
                          char filename_2[50], char filename_3[50])
{
 if ( (array_counters[0] != array_counters[1]) ||
      (array_counters[0] != array_counters[2]) ||
      (array_counters[1] != array_counters[2]) )
  {
   printf("%s\n","**ERROR: ARRAY COUNTS DISAGREE");
   printf("%s %s %s %d %s %s %s %d %s %s %s %d\n",
                 "**ERROR:",filename_1," Count=",array_counters[0]," - ",
                             filename_2," Count=",array_counters[1]," - ",
                             filename_3," Count=",array_counters[2]);
   printf("%s \n","**EXITING");
   exit(0);
  }
}



/* ===========================================================================
;
;  Name: check_files
;
;  Description: Check to see if output files already exist.
;               If they already exist, delete them.
;
;  Inputs: Names of output files.
;
;  Outputs: None.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   01/18/00   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>
#include <string.h>
#include <unistd.h>

/* ---------------------------------------------------------------------- */
void check_files(char output_file_names[4][50])
{
 int i;
 
 for (i = 0; i <= 3; ++i)
  {
    /* F_OK = Check existence of file */
   if(!access(output_file_names[i],F_OK))
    {
     printf("%s %s\n","Deleteing old: ",output_file_names[i]);
     unlink(output_file_names[i]);
    } 
  }
}





















/* ===========================================================================
;
;  Name: close_input_files
;
;  Description: close the input files.
;
;  Inputs: file pointers to input files
;
;  Outputs: None.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/15/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>

void close_input_files(FILE *input_files_ptr[3])
{
 int i;

 for (i = 0; i <= 2; ++i)
  fclose(input_files_ptr[i]);
}
               
/* ===========================================================================
;
;  Name: close_output_files
;
;  Description: close the output files.
;
;  Inputs: file pointers to output files
;
;  Outputs: None.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/15/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>

void close_output_files(FILE *output_files_ptr[4])
{
 int i;

 for (i = 0; i <= 3; ++i)
  fclose(output_files_ptr[i]);
}
               
/* ===========================================================================
;
;  Name: coord_xform
;
;  Description: Does transformation to various coord systems.
;
;  Inputs: b field and time data
;
;  Outputs: The b field transformed to the various coord. sys.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/15/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>
#include <math.h>

#include "mag_eros_products_data_types.h"
#include "mag_eros_products_prototypes.h"
#include "spice2.h"
#include "spice_id.h"

/* ---------------------------------------------------------------------- */
void coord_xform(int year, int doy, long int array_counter, 
                 struct time_data proc_time,
         struct time_errors *proc_time_errors,
                 struct b_data b_mrg,
                 struct b_conv_data *b_J2000,
                 struct position_data *eros_pos,
                 struct position_data *near_pos,
                 struct position_data *earth_pos,
                 struct velocity_data *eros_vel,
                 struct velocity_data *near_vel,
                 struct velocity_data *earth_vel,
         struct unit_vectors_J2000 *ABF_unit_vec_J2000,
                 struct b_conv_data *b_ABF,
                 struct position_data *dr_neros_ABF,
                 struct b_conv_data *b_SUN,
                 struct position_data *dr_neros_SUN)

{
 long int i;
 int j,m,n;
 double nmet,tcmat[3][3],cmat[3][3],av[3];
 int errst,handle;
 double ET,astermx[3][3],tastermx[3][3];
 double eros_state[6],near_state[6],earth_state[6];
 int yrday;
 int temp_e;
 double b_old[3],b_J2000_temp[3];
 double ex_ABF[3] = {1,0,0};
 double ez_ABF[3] = {0,0,1};
 double e1_sun[3], e2_sun[3], e3_sun[3],mage3_sun;
 double b_SUN_temp[3];
 double temp_calc;
 double ex_J2000[3], ez_J2000[3];
 double dr_neros_J2000[3]; /* position differences between near and eros in J2K*/
 double dr_neros_ABF_temp[3]; /* position differences between near and eros in ABF*/
 double e_veros[3], magv_eros;
 double v_eros[3], e_reros[3]; 
 double magr_eros, r_eros[3];
 double b_ABF_temp[3];
 double dr_neros_SUN_temp[3]; /* position differences between near and eros in SUN*/

 void mxv(double [][3],double[],double[]);

 int metet_err = 1; /* Assume success to start with these 4 error flags. */
 int near_state_err = 1;
 int earth_state_err = 1;
 int eros_state_err = 1;

 /* ------------------------------------------------------------ */

 for (i = 0; i <= array_counter; ++i)
  {
    /* J2000 calculations */

   metet_err = metet(NEAR, proc_time.MET[i], &ET, &errst);
   proc_time_errors->metet_error[i] = metet_err;
   if (metet_err) /* no error */
    {
     near_state_err = posvel(ET, NEAR, near_state, &errst);
     near_pos->error[i] = near_state_err;
     earth_state_err = posvel(ET, EARTH, earth_state, &errst);
     earth_pos->error[i] = earth_state_err;
     eros_state_err = posvel(ET, EROS, eros_state, &errst);
     eros_pos->error[i] = eros_state_err;
    }
   else
    {
     near_state_err = 0;
     earth_state_err = 0;
     eros_state_err = 0;
    }

   /* Divide distance by 1000 to report in 1000's of km.) */
   /* Velocity is reported in km/s. Multiply by 1000.0 to get m/s. */
  
   if (eros_state_err)
    {
     eros_pos->x[i] = eros_state[0] / 1000.0; 
     eros_pos->y[i] = eros_state[1] / 1000.0; 
     eros_pos->z[i] = eros_state[2] / 1000.0; 
     eros_vel->vx[i] = eros_state[3] * 1000.0; 
     eros_vel->vy[i] = eros_state[4] * 1000.0; 
     eros_vel->vz[i] = eros_state[5] * 1000.0; 
    }
   else
    {
     eros_pos->x[i] = -99.99; 
     eros_pos->y[i] = -99.99; 
     eros_pos->z[i] = -99.99; 
     eros_vel->vx[i] = -99.99; 
     eros_vel->vy[i] = -99.99; 
     eros_vel->vz[i] = -99.99; 
    }

   if (near_state_err)
    {
     near_pos->x[i] = near_state[0] / 1000.0; 
     near_pos->y[i] = near_state[1] / 1000.0; 
     near_pos->z[i] = near_state[2] / 1000.0; 
     near_vel->vx[i] = near_state[3] * 1000.0; 
     near_vel->vy[i] = near_state[4] * 1000.0; 
     near_vel->vz[i] = near_state[5] * 1000.0; 
    }
   else
    {
     near_pos->x[i] = -99.99; 
     near_pos->y[i] = -99.99; 
     near_pos->z[i] = -99.99; 
     near_vel->vx[i] = -99.99; 
     near_vel->vy[i] = -99.99; 
     near_vel->vz[i] = -99.99; 
    }

   if (earth_state_err)
    { 
     earth_pos->x[i] = earth_state[0] / 1000.0; 
     earth_pos->y[i] = earth_state[1] / 1000.0; 
     earth_pos->z[i] = earth_state[2] / 1000.0;       
     earth_vel->vx[i] = earth_state[3] * 1000.0; 
     earth_vel->vy[i] = earth_state[4] * 1000.0; 
     earth_vel->vz[i] = earth_state[5] * 1000.0;
    }
   else
    {
     earth_pos->x[i] = -99.99; 
     earth_pos->y[i] = -99.99; 
     earth_pos->z[i] = -99.99;
     earth_vel->vx[i] = -99.99; 
     earth_vel->vy[i] = -99.99; 
     earth_vel->vz[i] = -99.99;
    }

   /* Asteroid Body Fixed Coord. */
   if (near_state_err && eros_state_err)
    {
     for (j = 0; j < 3; j++)
      dr_neros_J2000[j] = near_state[j] - eros_state[j];   
    }
   else
    {
     for (j = 0; j < 3; j++)
      dr_neros_J2000[j] = -99.99;   
    }

   /* Sun coordinates; create unit vectores in J2000 */
  
   if (eros_state_err)
    {
     for (j = 3; j < 6; j++)
      v_eros[j-3] = eros_state[j];

     vmag(v_eros,&magv_eros);
 
     for (j = 0; j < 3; j++)
      e_veros[j] = v_eros[j] / magv_eros;

     for (j = 0; j < 3; j++)
      r_eros[j] = eros_state[j];

     vmag(r_eros,&magr_eros);
 
     for (j = 0; j < 3; j++)
      e1_sun[j] = -r_eros[j] / magr_eros;

     vcross(e_veros,e1_sun,e3_sun);
     vmag(e3_sun,&mage3_sun);

     for (j = 0; j < 3; j++)
      e3_sun[j] = e3_sun[j] / mage3_sun;

     vcross(e3_sun,e1_sun,e2_sun);  
    }   /* if (eros_state_err)  */

   nmet = metcmat2(NEAR_SC_BUS_PRIME, proc_time.MET[i], 1.0, cmat, av);

   if ( nmet < 0.0 )
    {
     proc_time_errors->cmat_error[i] = 0;
     printf("%s %lf %s %lf %s %d\n",
           "**ERROR: Trouble with: ",proc_time.MET[i]," nmet= ",nmet," i= ",i);
     b_J2000->bx[i] = -99.99;
     b_J2000->by[i] = -99.99;
     b_J2000->bz[i] = -99.99;
     b_J2000->b_tot[i] = -99.99;
     b_SUN->bx[i] = -99.99;
     b_SUN->by[i] = -99.99;
     b_SUN->bz[i] = -99.99;
     b_SUN->b_tot[i] - -99.99;
     dr_neros_SUN->x[i] = -99.99;
     dr_neros_SUN->y[i] = -99.99; 
     dr_neros_SUN->z[i] = -99.99;
     b_ABF->bx[i] = -99.99;
     b_ABF->by[i] = -99.99;
     b_ABF->bz[i] = -99.99;
     b_ABF->b_tot[i] = -99.99;
     ABF_unit_vec_J2000->ex1[i] = -99.99;
     ABF_unit_vec_J2000->ex2[i] = -99.99;
     ABF_unit_vec_J2000->ex3[i] = -99.99;
     ABF_unit_vec_J2000->ez1[i] = -99.99;
     ABF_unit_vec_J2000->ez2[i] = -99.99;
     ABF_unit_vec_J2000->ez3[i] = -99.99;
     dr_neros_ABF->x[i] = -99.99;
     dr_neros_ABF->y[i] = -99.99; 
     dr_neros_ABF->z[i] = -99.99;
    }
   else
    {
     proc_time_errors->cmat_error[i] = 1;
     b_old[0] = (-(1.0/sqrt(2.0))*b_mrg.bx[i]) - 
                ( (1.0/sqrt(2.0))*b_mrg.by[i]);
     b_old[1] = ( (1.0/sqrt(2.0))*b_mrg.bx[i]) - 
                ( (1.0/sqrt(2.0))*b_mrg.by[i]);
     b_old[2] = b_mrg.bz[i];

     /* Take the transpose of the c matrix provided by SPICE. */
     for (m = 0; m <= 2; ++m)
      for (n = 0; n <= 2; ++n) tcmat[m][n] = cmat[n][m];

     /* Convert the b field data tp J2000. */
     mxv(tcmat,b_old,b_J2000_temp);   

     b_J2000->bx[i] = b_J2000_temp[0];
     b_J2000->by[i] = b_J2000_temp[1];
     b_J2000->bz[i] = b_J2000_temp[2];
     b_J2000->b_tot[i] = sqrt(pow(b_J2000_temp[0],2.0) + 
                              pow(b_J2000_temp[1],2.0) +
                              pow(b_J2000_temp[2],2.0));
    
    /* Project J2000 field along e1,e2,e3 directions. */   

     if (eros_state_err)
      {
       vdot(b_J2000_temp,e1_sun,&temp_calc);
       b_SUN_temp[0] = temp_calc;
       vdot(b_J2000_temp,e2_sun,&temp_calc);
       b_SUN_temp[1] = temp_calc;
       vdot(b_J2000_temp,e3_sun,&temp_calc);
       b_SUN_temp[2] = temp_calc;

       b_SUN->bx[i] = b_SUN_temp[0];
       b_SUN->by[i] = b_SUN_temp[1];
       b_SUN->bz[i] = b_SUN_temp[2];
       b_SUN->b_tot[i] = sqrt(pow(b_SUN_temp[0],2.0) + 
                              pow(b_SUN_temp[1],2.0) +
                              pow(b_SUN_temp[2],2.0));

       vdot(dr_neros_J2000,e1_sun,&temp_calc);
       dr_neros_SUN_temp[0] = temp_calc;
       vdot(dr_neros_J2000,e2_sun,&temp_calc);
       dr_neros_SUN_temp[1] = temp_calc;
       vdot(dr_neros_J2000,e3_sun,&temp_calc);
       dr_neros_SUN_temp[2] = temp_calc;

       dr_neros_SUN->x[i] = dr_neros_SUN_temp[0];
       dr_neros_SUN->y[i] = dr_neros_SUN_temp[1]; 
       dr_neros_SUN->z[i] = dr_neros_SUN_temp[2];       
      }  /* if (eros_state_err) */
     else
      {
       b_SUN->bx[i] = -99.99;
       b_SUN->by[i] = -99.99;
       b_SUN->bz[i] = -99.99;
       b_SUN->b_tot[i] - -99.99;
       dr_neros_SUN->x[i] = -99.99;
       dr_neros_SUN->y[i] = -99.99; 
       dr_neros_SUN->z[i] = -99.99;
      }

     /*  Asteroid Body Fixed coordinates. */

     metet_err = metet(NEAR,proc_time.MET[i],&ET,&errst);
     proc_time_errors->metet_error[i] = metet_err;
     if (metet_err)
      {
       temp_e = EROS;
       bodmat(&temp_e,&ET,astermx);

       mxv(astermx,b_J2000_temp,b_ABF_temp);

       b_ABF->bx[i] = b_ABF_temp[0];
       b_ABF->by[i] = b_ABF_temp[1];
       b_ABF->bz[i] = b_ABF_temp[2];
       b_ABF->b_tot[i] = sqrt(pow(b_ABF_temp[0],2.0) + 
                              pow(b_ABF_temp[1],2.0) +
                              pow(b_ABF_temp[2],2.0));

     /* Asteroid Body Fixed unit vectors in J2000 */

     /* Take the transpose of the matrix. */
       for (m = 0; m <= 2; ++m)
        for (n = 0; n <= 2; ++n) tastermx[m][n] = astermx[n][m];

       mxv(tastermx,ex_ABF,ex_J2000);
       mxv(tastermx,ez_ABF,ez_J2000);
    
       ABF_unit_vec_J2000->ex1[i] = ex_J2000[0];
       ABF_unit_vec_J2000->ex2[i] = ex_J2000[1];
       ABF_unit_vec_J2000->ex3[i] = ex_J2000[2];
       ABF_unit_vec_J2000->ez1[i] = ez_J2000[0];
       ABF_unit_vec_J2000->ez2[i] = ez_J2000[1];
       ABF_unit_vec_J2000->ez3[i] = ez_J2000[2];

     /* Vector position from eros to near in ABF. */

       mxv(astermx,dr_neros_J2000,dr_neros_ABF_temp);   
    
       dr_neros_ABF->x[i] = dr_neros_ABF_temp[0];
       dr_neros_ABF->y[i] = dr_neros_ABF_temp[1]; 
       dr_neros_ABF->z[i] = dr_neros_ABF_temp[2];       
      }  /* if (metet_err)  */
     else
      {
       b_ABF->bx[i] = -99.99;
       b_ABF->by[i] = -99.99;
       b_ABF->bz[i] = -99.99;
       b_ABF->b_tot[i] = -99.99;
       ABF_unit_vec_J2000->ex1[i] = -99.99;
       ABF_unit_vec_J2000->ex2[i] = -99.99;
       ABF_unit_vec_J2000->ex3[i] = -99.99;
       ABF_unit_vec_J2000->ez1[i] = -99.99;
       ABF_unit_vec_J2000->ez2[i] = -99.99;
       ABF_unit_vec_J2000->ez3[i] = -99.99;
       dr_neros_ABF->x[i] = -99.99;
       dr_neros_ABF->y[i] = -99.99; 
       dr_neros_ABF->z[i] = -99.99;
      }

    }   /* if ( nmet < 0.0 ) */

 } /* for (i = 0; i <= array_counter; ++i) */

}

double dmin(double var1, double var2)
{
 double fpMin;

 if (var1 <= var2 )
  fpMin = var1;
 else
  fpMin = var2;
  
 return fpMin;
}
/* ===========================================================================
;
;  Name: get_input_file_list
;
;  Description: Get filename of file containing list of files to
;               process.
;
;  Inputs: None.
;
;  Outputs: Filename
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   01/30/97   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>
#include <string.h>

void  get_input_file_list(char input_file_list[75])
{
 printf("%s","Enter file containing list of files: ");
 scanf("%s",input_file_list);
}
/* ===========================================================================
;
;  Name: get_year_doy
;
;  Description: Grab year and doy from file name.
;
;  Inputs: Input file name.
;
;  Outputs: year and day of year
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   02/02/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>
#include <string.h>
#include <stdlib.h>

void  get_year_doy(char b_file_name[50], int *year, int *doy)
{
 char temp_year[3], temp_doy[4];
 int yy;

 strcpy(temp_year,"  ");
 strcpy(temp_doy,"   ");

 temp_year[0] = b_file_name[1];
 temp_year[1] = b_file_name[2];
  temp_doy[0] = b_file_name[3];
  temp_doy[1] = b_file_name[4];
  temp_doy[2] = b_file_name[5];

 yy = atoi(temp_year);

 if (yy == 0)
  yy = 2000;
 else
  yy = yy + 1900;

 *year = yy;
  *doy = atoi(temp_doy);
}

/* ==========================================================================
;
;  Name: apply_heater_correction.c
;
;  Description: Apply heater ON/OFF corrections to spacecraft mag. data.
;
;  Inputs: 
;
;  Outputs:
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/00     Original code
;
; $Id:$
;
; $Log: $
========================================================================== */

float by_htr(double t, double ton, double toff);

void heater_correction(double tb[],float bx[],float by[],float bz[],
                       double th0[],double th1[],int nb,int nh)
{
 int i,j;
 float dby;
 double ton, toff;

 j = 0;
 while( (tb[0] > th1[j]) && (j <= nh) )  ++j;

 while ( (i <= nb) && (j <= nh) ) 
  {
   i=0;
   while( (tb[i] < th0[j]) && (i <= nb) ) ++i;

   if ( i <= nb )
    {   
     ton = th0[j];
     toff = th1[j];

     while ( (tb[i] < th1[j]) && (i <= nb) )
      {
       dby = by_htr(tb[i],ton,toff);
       bx[i] = bx[i] - 0.1/5.7*dby;
       by[i] = by[i] - dby;
       bz[i] = bz[i] + 0.77/5.7*dby;
       ++i;
      } /* while ( (tb[i] < th1[j]) && (i <= nb) ) */   
    } /* if ( i <= nb ) */
   ++j;
  } /* while ( (i <= nb) && (j <= nh) ) */
}

/* ==========================================================================
;
;  Name: apply_heater_correction.c
;
;  Description: Apply heater ON/OFF corrections to spacecraft mag. data.
;
;  Inputs: 
;
;  Outputs:
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/00     Original code
;
; $Id:$
;
; $Log: $
========================================================================== */

float by_htr(double t, double ton, double toff);

void heater_correction(double tb[],float bx[],float by[],float bz[],
                       double th0[],double th1[],int nb,int nh)
{
 int i,j;
 int exit_proc = 0;
 int array_check1 = 1;
 int array_check2 = 1;
 int array_check3 = 1;
 float dby;
 double ton, toff;

 printf("%s\n","in heater_correction");

 j = 1;
 while (array_check3)
  {
   while (array_check1 && !exit_proc)
    {
     if (tb[1] > th1[j]) 
      ++j;
     else
      array_check1 = 0;
     if (array_check1)
      {
       if (j > nh)
        {
         exit_proc = 1; 
         array_check3 = 0;
        }
      }
    }

   if (!exit_proc)
    {
     i = 1;
     array_check1 = 1;
     while(array_check1 && !exit_proc)
      {
       if (tb[i] <= th0[j])
        ++i;    /* Increment tb if not within heater cycle. */
       else
        array_check1 = 0;
       if (array_check1)
    {
         if (i > nb)
      {
           exit_proc = 1;
           array_check3 = 0;
          }
        }
      }
    } /* if (!exit_proc) */

   if (!exit_proc)
    {
     ton = th0[j];
     toff = th1[j];

     array_check2 = 1;
     while (array_check2)
      {
       dby = by_htr(tb[i],ton,toff);
       bx[i] = bx[i] - 0.1/5.7*dby;
       by[i] = by[i] - dby;
       bz[i] = bz[i] + 0.77/5.7*dby;
       ++i;

       if (i > nb) 
        {
         exit_proc = 1;
         array_check2 = 0;
         array_check3 = 0;
        }
       else
        if (tb[i] > th1[j]) array_check2 = 0;
      }  /* while (array_check2) */

     if (!exit_proc)
      {
       ++j;
       if (j > nh) 
        {
         exit_proc = 1;
         array_check3 = 0;
        }
      }  /* if (!exit_proc) */
    }  /* if (!exit_proc) */
  }  /* while (array_check3) */
} 

/* ==========================================================================
;
;  Name: apply_heater_correction.c
;
;  Description: Apply heater ON/OFF corrections to spacecraft mag. data.
;
;  Inputs: 
;
;  Outputs:
;
;    Author          Date       Comment
;  ----------      --------   -----------
;  Bruce Toth       05/00     Original code
;
; $Id:$
;
; $Log: $
========================================================================== */

float by_htr(double t, double ton, double toff);

void heater_correction(double tb[],float bx[],float by[],float bz[],
                       double th0[],double th1[],int nb,int nh)
{
 int i,j;
 float dby;
 double ton, toff;

 printf("%s\n","in heater_correction");

 j = 0;
 while( (tb[0] > th1[j]) && (j < nh) )  ++j;
 if (j < nh) 
  {
   i=0;
   while( (tb[i] < th0[j]) && (i < nb) ) ++i;
   if (i < nb)
    {
     ton = th0[j];
     toff = th1[j];
     while ( (tb[i] < th1[j]) && (i < nb) )
      {
       dby = by_htr(tb[i],ton,toff);
       bx[i] = bx[i] - 0.1/5.7*dby;
       by[i] = by[i] - dby;
       bz[i] = bz[i] + 0.77/5.7*dby;
       ++i;
       }
     }
   }
 }

/* ===========================================================================
;
;  Name: mag_eros_products_main.c
;
;  Description: Produce the data files needed by the mag team for
;               the Eros encounter.
;
;  Inputs: byyddd.fix, bctpyyddd.mrg, and ctpyyddd.mrg
;
;  Outputs: b_SC_yyddd.dat, b_J2000.dat, b_ABF_yyddd.dat,b_ESO_yyddd.dat
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   09/15/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>
#include <string.h>
#include "spice2.h"

#include "mag_eros_products_data_types.h"
#include "mag_eros_products_prototypes.h"

/* ----------------------------------------------------------------------- */
main()
{
 char b_input_files[100][50];
 char ctp_input_files[100][50];
 char bctp_input_files[100][50];
 int num_b_files;
 char input_file_list[50];

/* ----------------------------------------------------------------------- */

 spiceinit("bruce.ker", -9);

/* Get the file containing the list if files to process from the user. */

 /*  get_input_file_list(input_file_list); */

  strcpy(input_file_list,"fix_files.lis");

/* Read the list of input files. */

 printf("%s \n","--> Reading list of input files.");

 read_input_file_list(input_file_list,&num_b_files,b_input_files);

/* Build the bctpyyddd.mrg and ctpyyddd.mrg file names based on the 
   byyddd.fix files read in. */

 printf("%s \n","--> Building bctp and ctp file names.");

 build_bctp_ctp_file_names(num_b_files,b_input_files,bctp_input_files,
                           ctp_input_files);

 process_b_bctp_ctp_files(num_b_files,b_input_files,bctp_input_files,
                          ctp_input_files);
}

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/*
Purpose:
Merges two time series with different and possibly erratic sampling intervals
choosing only those points which are within dt_max of each other from the two
series.  Only nearest points are paired.  Developed for use with NEAR MAG and
CTP data.  This task was first implemented in IDL using the 'where' statement. 
This  proved oppressively time consuming so a compiled version was required to
achieve  the necessary increase in execution speed.

Authors: B. J. Anderson and B. A. Toth, JHU/APL

Originated: October 19, 1999.

General description:
The two time series have sample times in seconds given by the arrays tb and
tc. 
The tb array is indexed using "i" while tc is indexed using "j".  The arrays 
are assumed to be monotonically increasing.  The routine first checks to see
which time  array has the larger next time step.  It then scans through the
more finely sampled time series to identify the two points that straddle the
morecoarsely sampled time series.  The closer of these two points is then
used to check for a match within dt_max seconds.  If the identified match is
a new time, it is added to the merged time and index lists.  Execution then
continues considering the next time steps.

Restrictions:
dt_max < 10,000.
tb, ti assumed to be monotonically increasing.
*/

/* ===================================================================== */
double dmin(double var1, double var2);

void merge_data(double dt_max,double tb[600000],long int na,
                double tc[600000],long int nb,int ri1[600000],
                int ri2[600000],long int *k) 
{
 long int i = 0;  /* Initialize indices and loop condition flags. */
 long int j = 0; 
 int flag = 0;
 int flag2 = 0;
 double da,db,d0,dm1;

 double tm[600000];
 int iii = 0;

 *k = 0;

 while (flag < 2)  /* flag ==2 is the exit condition when the end of either tb
              or tc is encountered  */
  {
   if (i < na)              /* da is next tb time step.  (Previous if this is 
                              the last point.) */
    da = tb[i+1] - tb[i];
   else
    da = tb[i] - tb[i-1];

   if (j < nb)              /* db is next tc time step.  (Previous if this is 
                               the last point.) */
    db = tc[j+1] - tc[j];
   else
    db = tc[j] - tc[j-1];

   if ( (da - db) > 0.001 )  /*  Identify the more finely sampled time series.

                 The > 0.001 condition ensures that tb is
                 treated is the more finely sampled time series
                 if the time steps are equal.  For NEAR this
                 was important because the CTP time series, tc,
                 occassionally had single point drop-outs
                 whereas the tb (MAG) time gaps are generally
                 at least 44 points (one packet) long.  It
                 therefore was cleaner to keep looping over tb
                 if the time steps were equal (happens for 1-s
                 CTP real time telemetry).

                 True implies tc is finer, so loop over j (tc
                 index). False implies tb is finer, so loop
                 over i (tb  index). */
    {
      flag = 0;
      while (flag == 0)      /*  tc-loop   */
       {
    if (j <= nb)
     {
      if (tc[j] >= tb[i])   /*  Increment j until tc[j] and tc[j-1] 
                                straddle tb[i]  */
       flag = 1;
      else
       ++j;
     } 
        else
     flag = 2;                /* Exit condition on end of tc array.  */
       }  /* while (flag == 0) */

      if (flag == 1)
       {
    d0 = fabs(tc[j]-tb[i]);   /* Construct abs. time differences between 
                                 tb[i] and bracketing tc points.  */
        if (j == 0)
     dm1 = 10000;             /* Default if first j point.  */
        else
     dm1 = fabs(tc[j-1]-tb[i]);

    if (dmin(d0,dm1) < dt_max)  /* Check if minimum time difference is less
                                   than dt_max.  */
     {
      if (*k > 0)
       {
        if (fabs(tb[i] - tm[*k-1]) > dmin(d0,dm1))
                    /*  Check to make sure current tb[i] is a new point. */
         {
           if (d0 < dm1)      /*  Assign to closer of j or j-1.  */
        {
         tm[*k] = tb[i];
         ri1[*k] = i;
             ri2[*k] = j;
                }
           else
        {
         tm[*k] = tb[i];
                 ri1[*k] = i;
             ri2[*k] = j-1;
                }
               ++*k;   /* Increment output merged array index after every new 
                      assignment.  */
         }
           }
      else
       {
        if (d0 < dm1)    /*  If first k, skip the new point check.  */
         {
          tm[*k] = tb[i];
          ri1[*k] = i;
          ri2[*k] = j;
             }
        else
         {
          tm[*k] = tb[i];
          ri1[*k] = i;
          ri2[*k] = j-1;
             }
        ++*k;
       } /* if (*k > 0) */
     } /* if (dmin(d0,dm1) < dt_max) */
    ++i;          /*  Increment to next i, tb index.  
                      (More coarsely sampled.)  */
        flag2 = 0;
        while(flag2 == 0)
     {
      if( (j+1) <= nb )
         /*  Now increment j up so that the jth point in tc is the last
             before tb[i]. Due to the vagaries of switching between the
             i and j loops, this guards against time gaps in the stream
             that had been treated as more finely sampled.  Although 
             not obvious, if there is a gap in j, the routine skips to
             the i-loop for one iteration and then returns to the
             j-loop.  On and subsequent to that return, the routine
             will  repeat entries with the same j index until the i
             index catches up to the gap in the tb time series.  Hence,
             this little loop.  */ 
       {
         if ( tb[i] <= tc[j+1] )
          flag2 = 1;
             else
          ++j;
           }
          else
       flag2 = 1;
         } /*  while(flag2 == 0) */
       } /* if (flag == 1) */
    } /* if (da > db) */
   else  
    {
     flag = 0;
     while (flag == 0 )      /*  tb-loop   */
      { 
       if (i <= na)
    {
     if (tc[j] <= tb[i])   /*  Increment i until tb[i] and tb[i-1] straddle
                               tc[j]  */
      flag = 1;
     else
      ++i;
    }
       else
    flag = 2;                /* Exit condition on end of tb array.  */
      } /* while(flag == 0) */
     if (flag == 1)
      {
    d0 = fabs(tc[j] - tb[i]);   /* Construct abs. time differences between 
                                   tc[j] and bracketing tb points.  */
        if (i == 0)
     dm1 = 10000;               /* Default if first i point. */
        else
         dm1 = fabs(tc[j] - tb[i-1]);

       if (dmin(d0,dm1) < dt_max)  /* Check if minimum time difference is less 
                                      than dt_max. */
    {
     if (*k > 0)
      {
       if (d0 < dm1)        /*  Assign to closer of i or i-1. */
        {
         if (fabs(tb[i] - tm[*k-1]) > dmin(d0,dm1))
               /*  Check to make sure current tb[i] is a new point. 
                   This loop structure is different than the j-loop
               because different tb's are involved and separate
               checks need to be performed for each case, i or i-1.*/
          {
           tm[*k] = tb[i];
           ri1[*k] = i;
           ri2[*k] = j;
               ++*k;       /* Increment output merged array index after every
                          new assignment.  */
              }
            }
       else
        {
         if (fabs(tb[i-1] - tm[*k-1]) > dmin(d0,dm1))
                /*  Check to make sure current tb[i-1] is a new point. */
          {
               tm[*k] = tb[i-1];
           ri1[*k] = i-1;
           ri2[*k] = j;
           ++*k;
              }
            }
          } /* if (*k > 0) */
     else
      {
       if (d0 < dm1)    /*  If first k, skip the new point check.  */
        { 
         tm[*k] = tb[i];
         ri1[*k] = i;
         ri2[*k] = j;
            }
       else
        {
         tm[*k] = tb[i-1];
         ri1[*k] = i-1;
         ri2[*k] = j;
            }
       ++*k;
          } /* if (*k > 0) */
    } /* if (dmin(d0,dm1) < dt_max) */
        ++j;                 /*  Increment to next j, i.e. tc index.
                             (More coarsely sampled.)  */
      } /* if (flag == 1) */
     flag2 = 0;
     while( flag2 == 0 )
      {
       if ( (i+1) <= na )
         /*  Now increment i up so that the ith point in tb is the last
             before tc[j]. The rationale for this check is the same as
             that given in  the j-loop.  */
    {
      if ( tc[j] <= tb[i+1] )
       flag2 = 1;
          else
       ++i;
        } 
       else
    flag2 = 1;
      } /* while( flag2 == 0 ) */
    } /* else of if (da > db) */
    if (i >= na) flag = 2;
    if (j >= nb) flag = 2;
  }  /* while (flag < 2) */
}
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double dmin(double var1, double var2);

void merge_data(double dt_max,double tb[600000],long int na,
                double tc[600000],long int nb,int ri1[600000],
                int ri2[600000],long int *k) 
{
 long int i = 0; 
 long int j = 0; 
 int flag = 0;
 int flag2 = 0;
 double da,db,d0,dm1;

 double tm[600000];
 int iii = 0;

 *k = 0;

 while (flag < 2)
  {
   if (i < na)
    da = tb[i+1] - tb[i];
   else
    da = tb[i] - tb[i-1];

   if (j < nb)
    db = tc[j+1] - tc[j];
   else
    db = tc[j] - tc[j-1];

   if ( (da - db) > 0.001 )
    {
      flag = 0;
      while (flag == 0)
       {
    if (j <= nb)
     {
      if (tc[j] >= tb[i])
       flag = 1;
      else
       ++j;
     } 
        else
     flag = 2;
       }  /* while (flag == 0) */

      if (flag == 1)
       {
    d0 = fabs(tc[j]-tb[i]);
        if (j == 0)
     dm1 = 10000;
        else
     dm1 = fabs(tc[j-1]-tb[i]);

    if (dmin(d0,dm1) < dt_max)
     {
      if (*k > 0)
       {
        if (fabs(tb[i] - tm[*k-1]) > dmin(d0,dm1))
         {
           if (d0 < dm1)
        {
         tm[*k] = tb[i];
         ri1[*k] = i;
             ri2[*k] = j;
                }
           else
        {
         tm[*k] = tb[i];
                 ri1[*k] = i;
             ri2[*k] = j-1;
                }
               ++*k;
         }
           }
      else
       {
        if (d0 < dm1)
         {
          tm[*k] = tb[i];
          ri1[*k] = i;
          ri2[*k] = j;
             }
        else
         {
          tm[*k] = tb[i];
          ri1[*k] = i;
          ri2[*k] = j-1;
             }
        ++*k;
       } /* if (*k > 0) */
     } /* if (dmin(d0,dm1) < dt_max) */
    ++i;
        flag2 = 0;
        while(flag2 == 0)
     {
      if( (j+1) <= nb )
       {
         if ( tb[i] <= tc[j+1] )
          flag2 = 1;
             else
          ++j;
           }
          else
       flag2 = 1;
         } /*  while(flag2 == 0) */
       } /* if (flag == 1) */
    } /* if (da > db) */
   else  
    {
     flag = 0;
     while (flag == 0 )
      { 
       if (i <= na)
    {
     if (tc[j] <= tb[i])
      flag = 1;
     else
      ++i;
    }
       else
    flag = 2;
      } /* while(flag == 0) */
     if (flag == 1)
      {
    d0 = fabs(tc[j] - tb[i]);
        if (i == 0)
     dm1 = 10000;
        else
         dm1 = fabs(tc[j] - tb[i-1]);

       if (dmin(d0,dm1) < dt_max)
    {
     if (*k > 0)
      {
       if (d0 < dm1)
        {
         if (fabs(tb[i] - tm[*k-1]) > dmin(d0,dm1))
          {
           tm[*k] = tb[i];
           ri1[*k] = i;
           ri2[*k] = j;
               ++*k;
              }
            }
       else
        {
         if (fabs(tb[i-1] - tm[*k-1]) > dmin(d0,dm1))
          {
               tm[*k] = tb[i-1];
           ri1[*k] = i-1;
           ri2[*k] = j;
           ++*k;
              }
            }
          } /* if (*k > 0) */
     else
      {
       if (d0 < dm1)
        { 
         tm[*k] = tb[i];
         ri1[*k] = i;
         ri2[*k] = j;
            }
       else
        {
         tm[*k] = tb[i-1];
         ri1[*k] = i-1;
         ri2[*k] = j;
            }
       ++*k;
          } /* if (*k > 0) */
    } /* if (dmin(d0,dm1) < dt_max) */
        ++j;    
      } /* if (flag == 1) */
     flag2 = 0;
     while( flag2 == 0 )
      {
       if ( (i+1) <= na )
    {
      if ( tc[j] <= tb[i+1] )
       flag2 = 1;
          else
       ++i;
        } 
       else
    flag2 = 1;
      } /* while( flag2 == 0 ) */
    } /* else of if (da > db) */
    if (i >= na) flag = 2;
    if (j >= nb) flag = 2;
  }  /* while (flag < 2) */
}
static char metcmat2_id[] = "$Id: metcmat2.c,v 1.4 1998/08/04 17:29:51 shyong Exp $";

/*
 * Get C-matrix from CK at a given MET.
 *
 * Input:
 *   double met      : Mission Elapse Time
 *   double tol      : Time tolerance in second
 *
 * Output:
 *   double cmat[][3]: C-matrix
 *   double av[]     : Angular Velocity Rate
 *
 * Return:
 *   Spacecraft clock time upon success, otherwise -1.
 * 
 * The C-Kernel will be loaded automatically.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "spice_id.h"

#define SCLKLEN    80
#define CKDIR "/project/near/SDC_archive/SpiceOps/ops/.ck"
#define FACTOR 1000

typedef int LOGICAL;

extern void reset(void);
extern LOGICAL failed(void);
extern void ckgpav(int *, double *, double *, char *, double[][3],
                   double[], double *, LOGICAL *, int);

double metcmat2(int inst, double met, double tol, double cmat[][3], double av[])
{
    double sclkdp;
    int sc;
    double clkout;
    char sclkch[SCLKLEN];
    char ref[6];
    LOGICAL fnd;
    int errst;
    char *strptr, *tmpptr;

    sc = inst / 1000;
    if (!metsclkdp(sc, met, &sclkdp, &errst))
        return -1;
    tol *= FACTOR;
    if (loadck2(met) < 0)
        return -1;

    strcpy(ref, "J2000");
    ckgpav(&inst, &sclkdp, &tol, ref, cmat, av, &clkout, &fnd, strlen(ref));
    if (fnd) {
        scdecd(&sc, &clkout, sclkch, SCLKLEN);
        if (failed()) {
            reset();
            return (0);
        }
        strptr = strchr(sclkch, '/');
        if (strptr == NULL)
            strptr = sclkch;
        else
            strptr++;
        tmpptr = strptr;
        while (isdigit(*tmpptr))
            tmpptr++;
        *tmpptr = '\0';
        return (atof(strptr) / FACTOR);                /* SUCCESS */
    } else
        return -1;
}

/*
 * Load CK for a given met
 * 
 * Input:
 *   double met  :  Mission Elapse Time (e.g., 43989400)
 *
 * Output: None
 *
 * Return:
 *   SPICE loaded file's handle or -1 when load failed.
 */
int loadck2(double met)
{
    static int  init=1;
    static int  p=0;
    static int  phandle, pyrday=0;
    static double smet=-86400;
    char ckfile[120];
    int handle, len;
    char dir[100], *tmpdir;
    int year, doy, yrday;
    int  errst;
    char utc[40];
    char str[20];

    void cklpf(char *, int *, int);
    void unldck2(int);

    if (met < smet || met >= smet + 86400) {
        metutc(NEAR, met, utc, &errst);
        utc[4] = '\0';
        utc[8] = '\0';
        year = atol(utc);
        doy = atol(utc+5);
        yrday = year * 1000 + doy;
        
        if (pyrday != 0)
            unldck2(phandle);

        if (NULL != (tmpdir = getenv("NEAR_CK")))
            strcpy(dir, tmpdir);
        else
            strcpy(dir, CKDIR);

        sprintf(ckfile, "%s/near%07d.bc", dir, yrday);
        fprintf(stderr, "loadck2 %s ...\n", ckfile);
        len = strlen(ckfile);
        cklpf(ckfile, &handle, len);
        if (failed()) {
           fprintf(stderr, "loadck2 %s fail\n", ckfile);
           reset();
           return -1;
        }
        phandle = handle;
        pyrday = yrday;
        sprintf(str, "%d-%03d//", year, doy);
        utcmet(NEAR, str, &smet, &errst);
    }
    return phandle;
}

/*
 * Unload a CK file
 *
 * input:
 *   int handle : Handle of CK file to be unloaded.
 */
void unldck2(int handle)
{
    void ckupf(int *);

    ckupf(&handle);
    if (failed())
        reset();
}
static char metquat_id[] = "$Id: metquat.c,v 1.1.1.1 1997/08/28 14:13:50 shyong Exp $";

/*
 * Get quaternion from CK at a given MET
 *
 * Input:
 *   double   met    : Mission Elapse Time
 *   double   tol    : Time tolerance in second
 *
 * Output:
 *   double   quat[] : Quaternion (q0, q1, q2, q3)
 *   double   av[]   : Angular Velocity Rate
 */

typedef int LOGICAL;

extern double metcmat2(int, double, double, double [][3], double []);
extern void m2q(double [][3], double []);

void metquat2(int inst, double met, double tol, double quat[], double av[])
{
    double cmat[3][3];
    double clkout;
    int    i;

    if (-1 == (clkout = metcmat2(inst, met, tol, cmat, av))) {
        for (i=0; i<4; i++)
            quat[i] = 0.0;
        for (i=0; i<3; i++)
            av[i] = 0.0;
    } else
        m2q(cmat, quat);
}
/* ===========================================================================
;
;  Name: open_input_files.c
;
;  Description: Open the requested byyddd.fix, bctpyyddd.mrg,a nd
;               ctpyyddd.mrg input files.
;
;  Inputs: Input file names.
;
;  Outputs: File ptr's and/or open error
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/15/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================== */

#include <stdio.h>

void open_input_files(char b_file[50],char bctp_file[50],char ctp_file[50],
                      FILE *input_files_ptr[2],int *open_error)
{
 int i;
 int all_files_ok = 1;
  
 for (i = 0; i <= 2; ++i)
  input_files_ptr[i] = NULL; /* no guarantee that the system will NULL it */
 
 input_files_ptr[0] = fopen(b_file,"r");
 input_files_ptr[1] = fopen(bctp_file,"r");
 input_files_ptr[2] = fopen(ctp_file,"r");

 if ( (input_files_ptr[0] == NULL) || (input_files_ptr[1] == NULL) ||
      (input_files_ptr[2] == NULL) )
  {
   all_files_ok = 0;
   *open_error = 1;
  }
 else
   *open_error = 0;

 /* If all files are not available, close files that have been opened
    because we will skip these files. */

 if (!all_files_ok)
  {
   for (i = 0; i <= 2; ++i)
    if ( input_files_ptr[i] != NULL) 
     fclose(input_files_ptr[i]);
  }

}
/* ===========================================================================
;
;  Name: open_output_files.c
;
;  Description: Open the output files.
;
;  Inputs: Input file names.
;
;  Outputs: File ptr's and/or open error
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/15/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================== */

#include <stdio.h>

void open_output_files(char output_file_names[4][50], 
                       FILE *output_file_ptrs[4], int *open_error)
{
 int i;
 int all_files_ok = 1;
  
 for (i = 0; i <= 3; ++i)
  {
   output_file_ptrs[i] = NULL; /* no guarantee that the system will NULL it */
   output_file_ptrs[i] = fopen(output_file_names[i],"a");
  }

 if ( (output_file_ptrs[0] == NULL) || (output_file_ptrs[1] == NULL) ||
      (output_file_ptrs[2] == NULL) || (output_file_ptrs[3] == NULL) )
  {
   all_files_ok = 0;
   *open_error = 1;
  }
 else
   *open_error = 0;

 /* If all files cannot be opened, close files that have been opened. */

 if (!all_files_ok)
  {
   for (i = 0; i <= 3; ++i)
    if ( output_file_ptrs[i] != NULL) 
     fclose(output_file_ptrs[i]);
  }
}
/* ===========================================================================
;
;  Name: print_input_file_open_errors
;
;  Description: Some files are not in the local directory, print
;               an error and exit.
;
;  Inputs: file pointers and file names
;
;  Outputs: None
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/15/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>

void print_input_file_open_errors(char b_file[50], char bctp_file[50],
                            char ctp_file[50])
{
 printf("%s %s %s %s %s %s\n",
  "**ERROR: Problem opening", b_file," and/or ",bctp_file," and/or ",
  ctp_file);
}
/* ===========================================================================
;
;  Name: print_output_file_open_errors
;
;  Description: Some files can't be opened in the local directory, print
;               an error and exit.
;
;  Inputs: file pointers and file names
;
;  Outputs: None
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   01/12/00   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>

void print_output_file_open_errors(char output_file_names[4][50])
{
 printf("%s %s %s %s %s %s %s %s\n",
  "**ERROR: Problem opening", output_file_names[0]," and/or ",
                              output_file_names[1]," and/or ",
                              output_file_names[2]," and/or ",
                              output_file_names[3]);
}
/* ===========================================================================
;
;  Name: process_b_bctp_ctp_files
;
;  Description: Process all the byyddd.fix, bctpyyddd.mrg, and
;               ctpyyddd.mrg files in order to produce the data
;               products required by the mag team for the Eros encounter.
;
;  Inputs: Number of files, byyddd.fix, bctpyyddd.mrg, and 
;          ctpyyddd.mrg file lists.
;
;  Outputs: b_SC_yyddd.dat, b_J2000_yyddd.dat, b_ABF_yyddd.dat,
;           and b_ESO_yyddd.dat
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/15/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>

#include "mag_eros_products_data_types.h"
#include "mag_eros_products_prototypes.h"

/* ---------------------------------------------------------------------- */
void process_b_bctp_ctp_files(int num_b_files,char b_input_files[100][50],
                              char bctp_input_files[100][50],
                              char ctp_input_files[100][50])
{
 int i;
 int open_input_files_error = 0;
 int open_output_files_error = 0;
 FILE *input_files_ptr[3];
 struct b_data b_mrg;
 struct bctp_data bctp_mrg; 
 struct ctp_data ctp_mrg; 
 struct time_data proc_time;
 struct time_errors proc_time_errors;

 struct b_conv_data b_J2000;
 struct position_data eros_pos;
 struct position_data near_pos;
 struct position_data earth_pos;
 struct velocity_data eros_vel;
 struct velocity_data near_vel;
 struct velocity_data earth_vel;
 struct unit_vectors_J2000 ABF_unit_vec_J2000;
 struct b_conv_data b_ABF;
 struct position_data dr_neros_ABF;
 struct b_conv_data b_SUN;
 struct position_data dr_neros_SUN;

 long int array_counters[3] = {0,0,0};
 int year, doy;
 char output_file_names[4][50];
 FILE *output_file_ptrs[4];
 int end_of_file = 0;
 int first_write = 1;

 for (i = 0; i <= num_b_files; ++i)
  {
    end_of_file = 0;
    open_input_files(b_input_files[i],bctp_input_files[i],ctp_input_files[i],
                     input_files_ptr,&open_input_files_error);
    
    if (!open_input_files_error)
     {
      printf("%s %s %s %s %s %s\n",
       "--> Processing: ",b_input_files[i]," and ",bctp_input_files[i],
       " and ",ctp_input_files[i]);
      /* Read the file headers - only once. */
      read_file_headers(input_files_ptr);

      /* Extract the year and doy from the file name. */
      get_year_doy(b_input_files[i],&year,&doy);

      /* Build the output file name. */
      printf("%s \n","--> Building output file names.");
      build_output_file_names(year,doy,output_file_names);

      /* Check if output files exist. If so, delete existing files. */
      printf("%s \n","--> Checking if output files exist. Will delete.");
      check_files(output_file_names);
     
      /* Open output files */
      open_output_files(output_file_names,output_file_ptrs,
                        &open_output_files_error);

      if (!open_output_files_error)
       {
    /* Process files in chunks based on array sizes. */
    while (!end_of_file)
     {
          /* Read the input files. */
          printf("%s \n","--> Reading input files.");
          read_input_files(input_files_ptr, &b_mrg, &bctp_mrg, &ctp_mrg, 
                           array_counters,&end_of_file);

          /* Check if b_data array counter and bctp_data array counter are the 
             same. If they disagree, flag an error and exit. */
          check_array_counters(array_counters, b_input_files[i],
                               bctp_input_files[i],ctp_input_files[i]);

          /* Calculate and load UTC and MET times. */
          printf("%s \n","--> Calculating UTC and MET times.");
          calc_load_times(array_counters[0],b_mrg,bctp_mrg,&proc_time,
                      &proc_time_errors);

          /* Here we transform the data to the various coord systems. */
          printf("%s \n","--> Transforming to various coord systems.");
 
          coord_xform(year,doy,array_counters[0],proc_time,&proc_time_errors,
                      b_mrg,&b_J2000,&eros_pos,&near_pos,&earth_pos,&eros_vel,
                      &near_vel,&earth_vel,&ABF_unit_vec_J2000,&b_ABF,
                      &dr_neros_ABF,&b_SUN,&dr_neros_SUN);

          /* Write the output file. */
          printf("%s \n","--> Writing the 4 output files.");
          write_output_files(output_file_ptrs,output_file_names,
                             array_counters[0],first_write,b_mrg,
                             bctp_mrg,ctp_mrg,proc_time,proc_time_errors,
                             b_J2000,
                             eros_pos,near_pos,earth_pos,eros_vel,near_vel,
                             earth_vel,ABF_unit_vec_J2000,b_ABF,dr_neros_ABF,
                             b_SUN,dr_neros_SUN);
          /* Set to FALSE so that the headers are not rewritten. */
          first_write = 0;
          /* Zero array counters to load next set of data. */
          zero_array_counters(array_counters);
         } /* while(!end_of_file) */
        close_input_files(input_files_ptr);
        close_output_files(output_file_ptrs);
       }
      else  /* (!open_output_files_error) */
       print_output_file_open_errors(output_file_names);
     }
    else
     print_input_file_open_errors(b_input_files[i],bctp_input_files[i],
                                  ctp_input_files[i]);
   
  }  /* for (i = 0; i <= num_b_files; ++i) */
}
/* ===========================================================================
;
;  Name: read_file_headers
;
;  Description: Read the file headers and "discard".
;
;  Inputs: file pointers
;
;  Outputs: None
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   01/12/00   Original code.
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>

void read_file_headers(FILE *input_files_ptr[3])
{
 int i;
 char temp_header[250];
  
 for (i = 0; i <= 2; i++)
  fscanf(input_files_ptr[i],"%[^\n]\n",temp_header);
}
/* ===========================================================================
;
;  Name: read_input_file_list
;
;  Description: Read the list of input files.
;
;  Inputs: Name of file containing list of input files.
;
;  Outputs: List of input files.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   01/30/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>

void read_input_file_list(char input_file_list[50],int *num_b_files,
                          char b_input_files[100][50])
{
 FILE *fp;
 int i = 0;

 fp = fopen(input_file_list,"r");

 if (fp != NULL)
  {
   while ( (fscanf(fp,"%s",b_input_files[i]) != EOF) && (i <= 99) ) ++i;
   --i;
   *num_b_files = i;
   fclose(fp);
  }
 else
  {
   printf("%s %s %s\n","**ERROR: TROUBLE OPENING ",input_file_list," EXITING");
   exit(0);
  }
}
/* ===========================================================================
;
;  Name: read_input_files
;
;  Description: Read in the input files byyddd.fix, bctpyyddd.mrg,
;               and ctpyyddd.mrg
;
;  Inputs: File pointers to input files.
;
;  Outputs: Arrays containing data from the input files.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/15/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>
#include <math.h>

#include "mag_eros_products_data_types.h"

void  read_input_files(FILE *input_files_ptr[3], struct b_data *b_mrg,
                       struct bctp_data *bctp_mrg, struct ctp_data *ctp_mrg,
                       long int array_counters[3], int *end_of_file)
{
 long int i,j;
 int t_day,t_hr,t_min;
 float t_sec,t_bx,t_by,t_bz;
 float t_diff;            /* dummy */
 int t_h1,t_h2,t_h3;      /* dummy */
 int t_ashnt1,t_dgshnt,t_slrarry;
 float factor1,factor2,factor3;
 float b_exp = 2.0;
 float temp1,temp1t,temp2t;
 int temp2,temp3,temp4,temp5,temp6,temp7,temp8;
 int continue_reading = 1;

 for (i = 0; i <= 2; ++i)
  {
   if ( i == 0 )     /* The byyddd.fix file. */
    {
     j = array_counters[i];
     continue_reading = 1;
     while (continue_reading)
      {
       if (j <= 49999)
    {
         fscanf(input_files_ptr[i],"%d%d%d%f%f%f%f%f%d",
            &t_day,&t_hr,&t_min,&t_sec,&temp1,
                &t_bx,&t_by,&t_bz,&temp2);
     *end_of_file = feof(input_files_ptr[i]);
     if (!*end_of_file)
      {
           b_mrg->day[j] = t_day;
           b_mrg->hr[j] = t_hr;
           b_mrg->min[j] = t_min;
           b_mrg->sec[j] = t_sec;
           b_mrg->bx[j] = t_bx;
           b_mrg->by[j] = t_by;
           b_mrg->bz[j] = t_bz;
           factor1 = pow(b_mrg->bx[j],b_exp);
           factor2 = pow(b_mrg->by[j],b_exp);
           factor3 = pow(b_mrg->bz[j],b_exp);
           b_mrg->b_tot[j] = sqrt(factor1 + factor2 + factor3);
       ++j;
          }   /* if (!*end_of_file) */
         else
      continue_reading = 0;
        }   /* if (j <= 49999) */
       else
    continue_reading = 0;
      }  /* while (continue_reading) */
 
     array_counters[i] = --j;
    }
   else if ( i == 1 )   /* The bctpyyddd.mrg file. */
    {
     j = array_counters[i];
     continue_reading = 1;
     while (continue_reading)
      {
       if (j <= 49999)
    {
         fscanf(input_files_ptr[i],"%d%d%d%f%f%f%f%f",
            &t_day,&t_hr,&t_min,&t_sec,&temp1,
                &t_bx,&t_by,&t_bz);
     *end_of_file = feof(input_files_ptr[i]);
     if (!*end_of_file)
      {
           bctp_mrg->day[j] = t_day;
           bctp_mrg->hr[j] = t_hr;
           bctp_mrg->min[j] = t_min;
           bctp_mrg->sec[j] = t_sec;
           bctp_mrg->bx[j] = t_bx;
           bctp_mrg->by[j] = t_by;
       ++j;
          }   /* if (!*end_of_file) */
         else
      continue_reading = 0;
        }   /* if (j <= 49999) */
       else
    continue_reading = 0;
      }  /* while (continue_reading) */
     array_counters[i] = --j;
    }
   else if ( i == 2 )  /* The ctpyyddd.mrg file. */
    {
     j = array_counters[i];
     continue_reading = 1;
     while (continue_reading)
      {
       if (j <= 49999)
    {
         fscanf(input_files_ptr[i],"%d%d%d%f%f%d%d%d%d%d%d%d%d%d%d%f%f",
            &temp2,&temp3,&temp4,&temp1,&t_diff,
                &t_h1,&t_h2,&t_h3,&t_ashnt1,&t_dgshnt,
                &t_slrarry,&temp5,&temp6,&temp7,&temp8,&temp1t,&temp2t);
     *end_of_file = feof(input_files_ptr[i]);
     if (!*end_of_file)
      {
           ctp_mrg->ASHNT1[j] = t_ashnt1;
           ctp_mrg->DGSHNT[j] = t_dgshnt;
           ctp_mrg->SLR_ARRY[j] = t_slrarry;
       ++j;
          }   /* if (!*end_of_file) */
         else
      continue_reading = 0;
        }   /* if (j <= 49999) */
       else
    continue_reading = 0;
      }  /* while (continue_reading) */
     array_counters[i] = --j;
    }
  } /* for i */
}


/* ===========================================================================
;
;  Name: vcross
;
;  Description: Calc cross product of 2 vectors.
;
;  Inputs: 2 vectors
;
;  Outputs: The cross product.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/15/98   Translated from Brian Anderson's FORTRAN code.
;
; $Id$
;
; $Log$
;
;========================================================================= */

void vcross(double v1[3],double v2[3], double v[3])
{
 v[0] = v1[1]*v2[2] - v1[2]*v2[1];
 v[1] = v1[2]*v2[0] - v1[0]*v2[2];
 v[2] = v1[0]*v2[1] - v1[1]*v2[0];
}
/* ===========================================================================
;
;  Name: vdot
;
;  Description: Calc dot product of 2 vectors.
;
;  Inputs: 2 vectors
;
;  Outputs: The dot product.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/15/98   Translated from Brian Anderson's FORTRAN code.
;
; $Id$
;
; $Log$
;
;========================================================================= */

void vdot(double v1[3],double v2[3], double *dot_prod)
{
 *dot_prod = v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
}
/* ===========================================================================
;
;  Name: vmag
;
;  Description: Calc magnitude of a vector.
;
;  Inputs: 1 vectors
;
;  Outputs: The magnitude of the vector.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/15/98   Translated from Brian Anderson's FORTRAN code.
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <math.h>

/* ----------------------------------------------------------------------- */
void vmag(double v1[3], double *v_mag)
{
 *v_mag = sqrt(v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2]);
}
/* ===========================================================================
;
;  Name: write_b_ABF_file
;
;  Description: Write the b_ABF_yyddd.dat output file.
;
;  Inputs: output file name, array counter, bctp & ctp data.
;
;  Outputs: None.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/17/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>

#include "mag_eros_products_data_types.h"

/*---------------------------------------------------------------------------*/
void write_b_ABF_file(FILE *output_file_ptr, long int array_counter,
                      int first_write,
                      struct time_data proc_time,struct b_data b_mrg,
                      struct b_conv_data b_ABF,
                      struct position_data dr_neros_ABF)
{
 long int i;

if (first_write)
 fprintf(output_file_ptr,"%s \n", 
 "UTC Year Day Hr Mn Sec MET Day Hr Mn Sec dMET BX BY BZ BTOT ABFx ABFy ABFz"); 

for (i = 0; i <= array_counter; ++i)
   fprintf(output_file_ptr,"%4d %4d %3d %3d %6.2f %4d %3d %3d %6.2f %7.2lf %9.2lf %9.2lf %9.2lf %9.2lf %15.4lf %15.4lf %15.4lf\n",
   proc_time.utc_year[i],proc_time.utc_doy[i],proc_time.utc_hr[i],
   proc_time.utc_min[i],proc_time.utc_sec[i],
   b_mrg.day[i],b_mrg.hr[i],b_mrg.min[i],b_mrg.sec[i],
   proc_time.MET_diff[i],b_ABF.bx[i],b_ABF.by[i], b_ABF.bz[i],
   b_ABF.b_tot[i],
   dr_neros_ABF.x[i],dr_neros_ABF.y[i],dr_neros_ABF.z[i]);
}
/* ===========================================================================
;
;  Name: write_b_ESO_file
;
;  Description: Write the b_ESO_yyddd.dat output file.
;               It was called  b_SUN_yyddd.dat output file.
;
;  Inputs: output file name, array counter, bctp & ctp data.
;
;  Outputs: None.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/17/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>

#include "mag_eros_products_data_types.h"

/*---------------------------------------------------------------------------*/
void write_b_ESO_file(FILE *output_file_ptr, long int array_counter,
                      int first_write,
                      struct time_data proc_time,struct b_data b_mrg,
                      struct b_conv_data b_SUN,
                      struct position_data dr_neros_SUN)
{
 long int i;

if (first_write)
 fprintf(output_file_ptr,"%s \n", 
 "UTC Year Day Hr Mn Sec MET Day Hr Mn Sec dMET BX BY BZ BTOT SUNx SUN SUNz"); 

for (i = 0; i <= array_counter; ++i)
   fprintf(output_file_ptr,"%4d %4d %3d %3d %6.2f %4d %3d %3d %6.2f %7.2lf %9.2lf %9.2lf %9.2lf %9.2lf %15.4lf %15.4lf %15.4lf\n",
   proc_time.utc_year[i],proc_time.utc_doy[i],proc_time.utc_hr[i],
   proc_time.utc_min[i],proc_time.utc_sec[i],
   b_mrg.day[i],b_mrg.hr[i],b_mrg.min[i],b_mrg.sec[i],
   proc_time.MET_diff[i],b_SUN.bx[i],b_SUN.by[i], b_SUN.bz[i],
   b_SUN.b_tot[i],
   dr_neros_SUN.x[i],dr_neros_SUN.y[i],dr_neros_SUN.z[i]);
}
/* ===========================================================================
;
;  Name: write_b_J2000_file
;
;  Description: Write the b_J2000_yyddd.dat output file.
;
;  Inputs: output file name, array counter, bctp & ctp data.
;
;  Outputs: None.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/17/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>

#include "mag_eros_products_data_types.h"

/*---------------------------------------------------------------------------*/
void write_b_J2000_file(FILE *output_file_ptr, long int array_counter,
                        int first_write,
                        struct time_data proc_time,
            struct time_errors proc_time_errors,
                        struct b_data b_mrg,
                        struct b_conv_data b_J2000,
                        struct position_data eros_pos,
                        struct position_data near_pos,
                        struct position_data earth_pos,
                        struct velocity_data eros_vel,
                        struct velocity_data near_vel,
                        struct velocity_data earth_vel,
                        struct unit_vectors_J2000 ABF_unit_vec_J2000)
{
 long int i; 

if (first_write)
 fprintf(output_file_ptr,"%s \n", 
"UTC Year Day Hr Mn Sec MET Day Hr Mn Sec dMET BX BY BZ BTOT R_Eros:x,y,z R_NEAR:x,y,z R_Earth:x,y,z V_Eros:vx,vy,vz V_NEAR:vx,vy,vz V_Earth:vx,vy,vz nx_Eros:nx1,nx2,nx3 nz_Eros:nz1,nz2,nz3 Error: metutc metet cmat state: eros near earth");

for (i = 0; i <= array_counter; ++i)
   fprintf(output_file_ptr,"%4d %4d %3d %3d %6.2f %4d %3d %3d %6.2f %7.2lf %9.2lf %9.2lf %9.2lf %9.2lf %16.7lf %16.7lf %16.7lf %16.7lf %16.7lf %16.7lf %16.7lf %16.7lf %16.7lf %10.2lf %10.2lf %10.2lf %10.2lf %10.2lf %10.2lf %10.2lf %10.2lf %10.2lf %10.5lf %10.5lf %10.5lf %10.5lf %10.5lf %10.5lf %d %d %d %d %d %d\n",
   proc_time.utc_year[i],proc_time.utc_doy[i],proc_time.utc_hr[i],
   proc_time.utc_min[i],proc_time.utc_sec[i],
   b_mrg.day[i],b_mrg.hr[i],b_mrg.min[i],b_mrg.sec[i],
   proc_time.MET_diff[i],b_J2000.bx[i],b_J2000.by[i], b_J2000.bz[i],
   b_J2000.b_tot[i],
   eros_pos.x[i],eros_pos.y[i],eros_pos.z[i],
   near_pos.x[i],near_pos.y[i],near_pos.z[i],
   earth_pos.x[i],earth_pos.y[i],earth_pos.z[i],
   eros_vel.vx[i],eros_vel.vy[i],eros_vel.vz[i],
   near_vel.vx[i], near_vel.vy[i], near_vel.vz[i],
   earth_vel.vx[i],earth_vel.vy[i],earth_vel.vz[i],
   ABF_unit_vec_J2000.ex1[i],ABF_unit_vec_J2000.ex2[i],
   ABF_unit_vec_J2000.ex3[i],ABF_unit_vec_J2000.ez1[i],
   ABF_unit_vec_J2000.ez2[i],ABF_unit_vec_J2000.ez3[i],
   proc_time_errors.metutc_error[i],proc_time_errors.metet_error[i],
   proc_time_errors.cmat_error[i],
   eros_pos.error[i],near_pos.error[i],earth_pos.error[i]);
}
/* ===========================================================================
;
;  Name: write_b_SC_file
;
;  Description: Write the b_SC_yyddd.dat output file.
;
;  Inputs: output file name, array counter, bctp & ctp data.
;
;  Outputs: None.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/17/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>

#include "mag_eros_products_data_types.h"

/*---------------------------------------------------------------------------*/
void write_b_SC_file(FILE *output_file_ptr, long int array_counter,
                     int first_write,
                     struct time_data proc_time, struct b_data b_mrg,
                     struct bctp_data bctp_mrg, struct ctp_data ctp_mrg)
{
 long int i;

if (first_write)
 fprintf(output_file_ptr,"%s \n", 
 "UTC Year Day Hr Mn Sec MET Day Hr Mn Sec dMET BX BY BZ BTOT ASHNT1 DGSHNT1 SLR_ARRAY BXCTP BYCTP BZCTP"); 

 for (i = 0; i <= array_counter; ++i)
   fprintf(output_file_ptr,"%4d %4d %3d %3d %6.2f %4d %3d %3d %6.2f %7.2lf %9.2f %9.2f %9.2f %9.2f %3d %3d %3d %7.2f %7.2f %7.2f\n",
   proc_time.utc_year[i], proc_time.utc_doy[i], proc_time.utc_hr[i],
   proc_time.utc_min[i], proc_time.utc_sec[i],
   b_mrg.day[i], b_mrg.hr[i], b_mrg.min[i], b_mrg.sec[i],
   proc_time.MET_diff[i],
   b_mrg.bx[i], b_mrg.by[i], b_mrg.bz[i], b_mrg.b_tot[i],
   ctp_mrg.ASHNT1[i], ctp_mrg.DGSHNT[i], ctp_mrg.SLR_ARRY[i],
   bctp_mrg.bx[i], bctp_mrg.by[i], bctp_mrg.bz[i]);
}
/* ===========================================================================
;
;  Name: write_output_files
;
;  Description: Write data to output files.
;
;  Inputs: output file names, array counter, b,bctp, & ctp data.
;
;  Outputs: None.
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   10/17/98   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */

#include <stdio.h>
#include <string.h>

#include "mag_eros_products_data_types.h"
#include "mag_eros_products_prototypes.h"

/* ----------------------------------------------------------------------- */
void write_output_files(FILE *output_file_ptrs[4],
                        char output_file_names[4][50],
                        long int array_counter,
                        int first_write,
                        struct b_data b_mrg,struct bctp_data bctp_mrg,
                        struct ctp_data ctp_mrg,struct time_data proc_time,
            struct time_errors proc_time_errors,
                        struct b_conv_data b_J2000,
                        struct position_data eros_pos,
                        struct position_data near_pos,
                        struct position_data earth_pos,
                        struct velocity_data eros_vel,
                        struct velocity_data near_vel,
                        struct velocity_data earth_vel,
                        struct unit_vectors_J2000 ABF_unit_vec_J2000,
                        struct b_conv_data b_ABF,
                        struct position_data dr_neros_ABF,
                        struct b_conv_data b_SUN,
                        struct position_data dr_neros_SUN)
{
   /* Write b_SC_yyddd.dat file. */
   printf("%s %s\n","Writing: ",output_file_names[0]);
   write_b_SC_file(output_file_ptrs[0],array_counter,first_write,proc_time,
                   b_mrg,bctp_mrg,ctp_mrg);

   /* Write b_J2000_yyddd.dat file. */

   printf("%s %s\n","Writing: ",output_file_names[1]);
   write_b_J2000_file(output_file_ptrs[1],array_counter,first_write,proc_time,
              proc_time_errors,b_mrg,
                      b_J2000,eros_pos,near_pos,earth_pos,eros_vel,near_vel,
                      earth_vel,ABF_unit_vec_J2000);
  
   /* Write b_ABF_yyddd.dat file. */

   printf("%s %s\n","Writing: ",output_file_names[2]);
   write_b_ABF_file(output_file_ptrs[2],array_counter,first_write,proc_time,
                    b_mrg,b_ABF,dr_neros_ABF);
 
   /* Write b_ESO_yydd.dat file. (was b_SUN) */

   printf("%s %s\n","Writing: ",output_file_names[3]);
   write_b_ESO_file(output_file_ptrs[3],array_counter,first_write,proc_time,
                    b_mrg,b_SUN,dr_neros_SUN);
}

/* ===========================================================================
;
;  Name: zero_array_counters
;
;  Description: zero out the array counters array
;
;  Inputs: array of array counters
;
;  Outputs: zeroed out array of array counters
;
;    Author       Date       Comment
;  ----------   --------   -----------
;  Bruce Toth   01/18/00   Original code
;
; $Id$
;
; $Log$
;
;========================================================================= */
void zero_array_counters(long int array_counters[3])
{
  int i;

  for (i = 0; i <= 2; ++i)
   array_counters[i] = 0;
}
struct b_data
 {
  int year[50000];
  int day[50000];
  int hr[50000];
  int min[50000];
  float sec[50000];
  float bx[50000];
  float by[50000];
  float bz[50000];
  float b_tot[50000];
 };

struct bctp_data
 {                
  int day[50000];
  int hr[50000];
  int min[50000];
  float sec[50000];
  float bx[50000];
  float by[50000];
  float bz[50000];
  float b_tot[50000];
 };

struct ctp_data
 {                
  int ASHNT1[50000];
  int DGSHNT[50000];
  int SLR_ARRY[50000];
 };

struct time_data
 {
  int utc_year[50000];
  int utc_doy[50000];
  int utc_hr[50000];
  int utc_min[50000];
  float utc_sec[50000];
  double MET[50000];
  double MET_diff[50000];
 };

struct time_errors
 {
  int metutc_error[50000];
  int metet_error[50000];
  int cmat_error[50000];
 };

struct position_data
 {
  double x[50000];
  double y[50000];
  double z[50000];
  int error[50000];    /* record state vector error */
 };

struct velocity_data
 {
  double vx[50000];
  double vy[50000];
  double vz[50000];
 };

struct unit_vectors_J2000
 {
  double ex1[50000];
  double ex2[50000];
  double ex3[50000];
  double ez1[50000];
  double ez2[50000];
  double ez3[50000];
 };

struct b_conv_data      /* Data converted to given coord sys. */ 
 {
  float bx[50000];
  float by[50000];
  float bz[50000];
  float b_tot[50000];
 };



void get_input_file_list(char input_file_list[50]);
void read_input_file_list(char input_file_list[50],int *num_b_files,
                          char b_input_files[100][50]);
void build_bctp_ctp_file_names(int num_b_files,char b_input_files[100][50],
                               char bctp_input_files[100][50],
                               char ctp_input_files[100][50]);
void process_b_bctp_ctp_files(int num_b_files,char b_input_files[100][50],
                              char bctp_input_files[100][50],
                              char ctp_input_files[100][50]);
void open_input_files(char b_file[50],char bctp_file[50],char ctp_file[50],
                      FILE *input_files_ptr[3],int *open_error);
void print_input_file_open_errors(char b_file[50], char bctp_file[50],
                                 char ctp_file[50]);
void print_output_file_open_errors(char output_file_names[4][50]);
void read_input_files(FILE *input_files_ptr[3], struct b_data *b_mrg,
                      struct bctp_data *bctp_mrg, struct ctp_data *ctp_mrg,
                      long int array_counters[3], int *end_of_file);
void close_input_files(FILE *input_files_ptr[3]);
void read_file_headers(FILE *input_files_ptr[3]);
void check_array_counters(long int array_counters[3], char filename_1[50], 
                          char filename_2[50], char filemame_3[50]);
void zero_array_counters(long int array_counters[3]);
void calc_load_times(long int array_count, struct b_data b_mrg,
                     struct bctp_data bctp_mrg, struct time_data *proc_time,
                     struct time_errors *proc_time_errors);
void calc_MET_time(int dday, int hhr, int mmin, double temp_sec,double *cbMET);
void get_year_doy(char b_file_name[50], int *year, int *doy);
void coord_xform(int year, int doy, long int array_counter, 
                 struct time_data proc_time,
         struct time_errors *proc_time_errors,
                 struct b_data b_mrg,
                 struct b_conv_data *b_J2000,
                 struct position_data *eros_pos,
                 struct position_data *near_pos,
                 struct position_data *earth_pos,
                 struct velocity_data *eros_vel,
                 struct velocity_data *near_vel,
                 struct velocity_data *earth_vel,
         struct unit_vectors_J2000 *ABF_unit_vec_J2000,
                 struct b_conv_data *b_ABF,
                 struct position_data *dr_neros_ABF,
                 struct b_conv_data *b_SUN,
                 struct position_data *dr_neros_SUN);
void calc_utc_hr_min_sec(double sod,int *utc_hr, int *utc_min, float *utc_sec);
void build_output_file_names(int year, int doy, char output_file_name[4][50]);
void check_files(char output_file_names[4][50]);
void open_output_files(char output_file_names[4][50],
                       FILE *output_file_ptrs[4],int *);
void write_output_files(FILE *output_file_ptrs[4],
                        char output_file_names[4][50],
                        long int array_counters,
            int first_write,
                        struct b_data b_mrg,struct bctp_data bctp_mrg,
                        struct ctp_data ctp_mrg,struct time_data proc_time,
            struct time_errors proc_time_errors,
                        struct b_conv_data b_J2000,
                        struct position_data eros_pos,
                        struct position_data near_pos,
                        struct position_data earth_pos,
                        struct velocity_data eros_vel,
                        struct velocity_data near_vel,
                        struct velocity_data earth_vel,
                        struct unit_vectors_J2000 ABF_unit_vec_J2000,
                        struct b_conv_data b_ABF,
                        struct position_data dr_neros_ABF,
                        struct b_conv_data b_SUN,
                        struct position_data dr_neros_SUN);
void write_b_SC_file(FILE *output_file_ptr, long int array_counter,
                     int first_write,
                     struct time_data proc_time, struct b_data b_mrg,
                     struct bctp_data bctp_mrg, struct ctp_data ctp_mrg);
void write_b_J2000_file(FILE *output_file_ptr, long int array_counter,
                        int first_write,
                        struct time_data proc_time, 
            struct time_errors proc_time_errors,
                        struct b_data b_mrg,
                        struct b_conv_data b_J2000,
                        struct position_data eros_pos,
                        struct position_data near_pos,
                        struct position_data earth_pos,
                        struct velocity_data eros_vel,
                        struct velocity_data near_vel,
                        struct velocity_data earth_vel,
                        struct unit_vectors_J2000 ABF_unit_vec_J2000);
void write_b_ABF_file(FILE *output_file_ptr, long int array_counter,
                      int first_write,
                      struct time_data proc_time,struct b_data b_mrg,
                      struct b_conv_data b_ABF,
                      struct position_data dr_neros_ABF);
void write_b_SUN_file(FILE *output_file_ptr, long int array_counter,
                      int first_write,
                      struct time_data proc_time,struct b_data b_mrg,
                      struct b_conv_data b_SUN,
                      struct position_data dr_neros_SUN);
void close_output_files(FILE *output_files_ptr[4]);
void vcross(double v1[3],double v2[3], double v[3]);
void vdot(double v1[3],double v2[3], double *dot);
void vmag(double v1[3], double *v_mag);
CC = cc -Aa +z -c

OBJS = merge_data.o dmin.o Callmerge_data.o heater_correction.o\
    by_htr.o a0fit.o a1fit.o a2fit.o a3fit.o Callheater_correction.o

pt_hp.so : $(OBJS)
    ld -b -o idlmerge_hp.sl $(OBJS)

merge_data.o : merge_data.c
    $(CC) -o merge_data.o merge_data.c -lm

dmin.o : dmin.c
    $(CC) -o dmin.o dmin.c

Callmerge_data.o : Callmerge_data.c
    $(CC) -o Callmerge_data.o Callmerge_data.c

heater_correction.o : heater_correction.c
    $(CC) -o heater_correction.o heater_correction.c

by_htr.o : by_htr.c
    $(CC) -o by_htr.o by_htr.c -lm

a0fit.o : a0fit.c
    $(CC) -o a0fit.o a0fit.c -lm

a1fit.o : a1fit.c
    $(CC) -o a1fit.o a1fit.c -lm

a2fit.o : a2fit.c
    $(CC) -o a2fit.o a2fit.c -lm

a3fit.o : a3fit.c
    $(CC) -o a3fit.o a3fit.c -lm

Callheater_correction.o : Callheater_correction.c
    $(CC) -o Callheater_correction.o Callheater_correction.c




CC = cc
CFLAGS = -Aa -D_HPUX_SOURCE -c

SDIR = /project/near/SDC_sw/PROD/sdc
SINC = -I$(SDIR)/include
SLIB = -lm -L$(SDIR)/lib -lspice2
SPDIR = /project/near/spice/naif/toolkit/lib
SPLIB = -L$(SPDIR)/ -lnear -lsupport -ltoolkit -lcl

.SUFFIXES:
.SUFFIXES: .c .o

.c.o:
    $(CC) $(CFLAGS) $(SINC) $<

OBJS = mag_eros_products_main.o \
open_input_files.o print_input_file_open_errors.o \
print_output_file_open_errors.o read_input_files.o \
read_input_file_list.o close_input_files.o read_file_headers.o\
build_bctp_ctp_file_names.o process_b_bctp_ctp_files.o \
check_array_counters.o calc_load_times.o calc_MET_time.o \
coord_xform.o get_year_doy.o build_output_file_names.o \
open_output_files.o write_output_files.o write_b_SC_file.o \
write_b_J2000_file.o write_b_ABF_file.o write_b_ESO_file.o \
close_output_files.o vcross.o vdot.o vmag.o check_files.o \
zero_array_counters.o metcmat2.o metquat2.o

mag_eros_products : $(OBJS)
    $(CC) -o $@ $(OBJS) $< $(SLIB) $(SPLIB)
--------------------------
.run magview_merge.pro
exit
---------------------------
     program mrg_corrections_v8
      
c
c
c  Filenames: bYYDDD.mrg       10 chars, format('b',a5,'.mrg')
c             bctpYYDDD.mrg    13 chars, format('bctp',a5,'.mrg')
c             ctpYYDDD.mrg     12 chars, format('ctp',a5,'.mrg')
c
c  Heritage: mrg_proc01.f  (For getting averages over days)
c
c  Pupose: to apply S/C corrections:
c    1.  Subtract BCTP
c    2.  Subtract ASHNT1 field.
c    3.  Subtract digital shunt steps.
c    4.  Subtract constant offset.
c    5.  Subtract range offsets for days before 578 (upload).
c    6.  Subtract solar array current contribution (1/22/98)
c    7.  Sensor orientation adjustment (2/20/98)
c
c  Digital shunt steps are subtracted relative to step #18.
c  Constant offset is evaluated for shunt step #18.
c  Range offsets are for instrument range 3.
c
c  Modified: 12/14-22/98
c     1.  New analog shunt model - six segment linear piecewise function
c     2.  Modified digital shunt steps.

c  Modified: 12/22/98
c     1.  No digital shunt steps - now differences in slopes
c         of arry current dependence.  (Modeling indicates shunt steps should be
c         invisible.)
c     3.  Incorporated better digital shunt current levels.
c

c  Modified: 1/4/98
c     1.  Scrub differences in slopes.  Now use shunt steps proportional to
c         total solar array current.
c
c
c  Modified: 2/13/00
c     1.  New digital shunt step coefficients derived and inserted.
c     2.  Offset proportional to electronics temperature correction added.
c         (Required access to bhsk files directly.)
c
c  Modified: 7/5/00
c     1.  Revised digital shunt step coefficients including data from >=100km
c         orbits.  Made significant difference in by coefs for steps 8 and 9.
c     2.  Added spacecraft offset correction depending on met_day from a 
c         lookup table.
c         Offsets determined by B.dB=0 method on 10 minute averages (despiked)
c         using a 27-day time span.
c
c  Modified: 3/19/01
c     1.  Added thruster valve correction to remove step from valve opening
c         that occurred during descent.  Step calibrated using 6 prior valve
c         operations.
c

      character*1 junk
      character*80 filelist
      character*5 yd
      character*10 bfile,outbfile
      character*12 ctpfile,scfld_file
      character*13 bctpfile,bhskfile
      real bsec,bt,bx,by,bz,bcsec,bct,bcx,bcy,bcz,csec,ct,
     &     dig_rat,dg_lowr(20),dg_uppr(20),
c     &     bx_ashnt(20),by_ashnt(20),bz_ashnt(20),
     &     bx_dshunt(20),by_dshunt(20),bz_dshunt(20),
     &     ox(8),oy(8),oz(8),
     &     arry_bx(20),arry_by(20),arry_bz(20),
     &     bxoffset(2),byoffset(2),bzoffset(2),
     &     thxy,thxz,cos_xy,sin_xy,cos_xz,sin_xz,
     &     bxtmp,bytmp,bztmp,
     &     xash,yash,zash,bxSC,bySC,bzSC,
     &     bash(3),boffset(3),barry(3),bdshunt(3),btelec(3),boff_SC(3)
      integer bday,bhr,bmn,rng,bcday,bchr,bcmn,cday,chr,cmn,
     &        inst_htr,sc_htr,prp_htr,ashnt1,dgshnt,slr_arry,
     &        id,scfld_unit,str_cam,nm_dpu,xg_dpu,aiu_1,
     &        dy_hsk,hr_hsk,mn_hsk,bhsk_unit,n_hsk,i,isc0,noff,isc1,k_sc0
      real tube_temp,probe_temp,
     &     brate,spikes,cal_on,toffset,rrange,s_hsk,dt_hsk,
     &     rman_rng,idcdc,tprobe,telect(14000),tdpu,t_interp,
     &     bx0_sc(500),by0_sc(500),bz0_sc(500),wt0,wt1
     &     bxTV,byTV,bzTV
      real*8 tmrg,tbhsk(14000),ttmp,t0_SC(500),xdy0_SC(500)

c  Versions 1&2.
c     
c      data dg_lowr/0.000,0.050,0.100,0.150,0.200,0.250,0.300,   ! guesss
c     &              0.350,0.394,0.439,0.485,0.529,0.570,0.616,   ! calibration values
c     &              0.668,0.721,0.774,0.826,0.878,               ! calibration values
c     &              0.950/                                       ! guess
c     
c      data dg_uppr/0.050,0.100,0.150,0.200,0.250,0.300,         ! guess
c     &              0.350,0.394,0.439,0.485,0.529,0.570,0.616,   ! calibration values
c     &              0.668,0.721,0.774,0.826,0.878,               ! calibration values
c     &              0.950,1.000/                                 ! guess
     
c  Version 3.          
      data dg_lowr/0.0000, 0.0495, 0.0990, 0.1484, 0.1979, 0.2474,
     &             0.2969, 0.3464, 0.3940, 0.4390, 0.4850, 0.5290,
     &             0.5700, 0.6160, 0.6680, 0.7210, 0.7740, 0.8260,
     &             0.8780, 0.9290/

      data dg_uppr/0.0495, 0.0990, 0.1484, 0.1979, 0.2474, 0.2969,
     &             0.3464, 0.3940, 0.4390, 0.4850, 0.5290, 0.5700,
     &             0.6160, 0.6680, 0.7210, 0.7740, 0.8260, 0.8780,
     &             0.9290, 1.0000/


c
c  Version 1 Analog shunt correction.  Obsolete.  Replaced with Version 2.
c  The real functions xash,yash,zash implement the Version 2 six element piecewise
c  linear fits that are used in Version 2.
c
c     
c      data bx_ashnt/7*0.0654,                                   ! avg. of cal. values
c     &                0.0530, 0.0641, 0.0639, 0.0634, 0.0757,    ! calibration values
c     &                0.0722, 0.0662, 0.0645,                    ! calibration values
c     &                5*0.0654/                                  ! avg. of cal. values
c     
c      data by_ashnt/7*0.0635,                                   ! avg. of cal. values
c     &                0.0921, 0.0766, 0.0693, 0.0703, 0.0587,    ! calibration values
c     &                0.0485, 0.0531, 0.0393,                    ! calibration values
c     *                5*0.0635/                                  ! avg. of cal. values
c     
c      data bz_ashnt/7*-0.0021,                                  ! avg. of cal. values
c     &                0.0044, 0.0174, 0.0096, 0.0030,-0.0231,    ! calibration values
c     &               -0.0180, 0.0071,-0.0175,                    ! calibration values
c     &                5*-0.0021/                                 ! avg. of cal. values
c     

c  Version 1: 
c      data bxoffset,byoffset,bzoffset/
c     &      86.6,  84.4,
c     &     -89.7, -69.2,
c     &     151.9, 160.5/
c         
c      data bx_dshunt,by_dshunt,bz_dshunt/
c     &     7*-4.67, -4.67, -6.02, -5.90, -5.57, -6.50,
c     &       -5.40, -4.24, -2.71, -5.32, -2.03, 3*0.0,
c     &     7*-3.66, -3.66, -1.94, -0.12,  0.70, -0.45,
c     &       -2.03,  2.15,  3.45, -1.74, -5.11, 3*0.0,
c     &     7*-1.91, -1.91, -3.07, -4.23, -4.91, -2.43,
c     &       -1.56, -4.62, -3.63, -1.66,  1.43, 3*0.0/

c  Version 2: 12/16/98
c      data bxoffset,byoffset,bzoffset/
c     &      81.9,  79.7,
c     &     -93.4, -72.9,
c     &     150.0, 158.6/
c         
c      data bx_dshunt,by_dshunt,bz_dshunt/
c     &     8*0.0, 0.21, 0.44, -0.07, 0.41, 3*0.33, 4.03, 4*6.13,
c     &     8*0.0, 2.01, 2.79,  1.14,-0.58, 3*1.85,-1.36, 4*2.72,
c     &     8*0.0,-0.94,-1.54, -0.18, 0.92, 3*1.58, 4.03, 4*2.77/


c  Version 3: 12/22/98
c      data bxoffset,byoffset,bzoffset/
c     &      81.9,  79.7,
c     &     -88.7, -68.2,
c     &     146.9, 155.5/
c         
c      data bx_dshunt,by_dshunt,bz_dshunt/
c     &     60*0.0/
c     
c      data arry_bx,arry_by,arry_bz/
c     & 8*-0.0198, -0.004965, -0.016297, -0.008384, -0.020320,
c     &  -0.025339, -0.019034, -0.025775, -0.040568, -0.013899,
c     &  -0.019098, 2*-0.0198,
c     &   8*0.1633,  0.155459,  0.154189,  0.169822,  0.155668,
c     &   0.139088,  0.195068,  0.201969,  0.170787,  0.143272,
c     &   0.169308,  2*0.1633,
c     &  8*-0.2594, -0.241036, -0.251677, -0.284205, -0.276618,
c     &  -0.246582, -0.275155, -0.282360, -0.257410, -0.243238,
c     &  -0.265082, 2*-0.2594/
     
c  Version 4: 12/22/98
      data bxoffset,byoffset,bzoffset/
     &      81.9,  79.7,
     &     -88.7, -68.2,
     &     146.9, 155.5/
        
c      data bx_dshunt,by_dshunt,bz_dshunt/
c     &     7*0.0,-0.0037,-0.0018,0.0018,0.0036,-0.0026,0.0043,
c     &     0.0039,0.0039,0.0039,0.0289,0.0395,0.0395,0.0395,
c     &     7*0.0,-0.0372,-0.0472,-0.0154,-0.0039,-0.0251,-0.0393,
c     &     0.013,0.013,0.013,-0.0084,0.0174,0.0174,0.0174,
c     &     7*0.0,0.0267,0.0321,0.0226,0.0132,0.0351,0.0533,
c     &     0.0325,0.0325,0.0325,0.0517,0.0447,0.0447,0.0447/
     
      data arry_bx,arry_by,arry_bz/
     &  20*-0.020,20*0.163,20*-0.259/
   
      data ox/0.0,0.0,0.0,0.0,-2.415,-1.565,-46.515,-34.015/
      data oy/0.0,0.0,0.0,0.0, 1.635, 1.640, 30.105, 15.055/
      data oz/0.0,0.0,0.0,0.0,-3.275,-2.135,-54.990,-35.630/
      
c      data arry_bx,arry_by,arry_bz/ -0.058,  0.158, -0.286/   ! Fit values to 97b period

      data thxy,thxz/1.76,0.22/   !  Sensor orientation correction angles.


c  Version 5: 02/12/00  Digital shunt step revision only.


      data bx_dshunt,by_dshunt,bz_dshunt/
     &     8*-0.01338,-0.00823,-0.00381,-0.00014,-0.00846,0.0,
     &     0.00323,0.01260,-0.01155,0.01157,3*0.02280,
     &     8*0.03105,0.00887,0.03742,0.04927,0.02340,0.0,
     &     0.06126,0.06719,0.03205,0.00575,3*0.03260,
     &     8*-0.04375,-0.02438,-0.03195,-0.04400,-0.02305,0.0,
     &     -0.01983,-0.03545,-0.00768,0.01908,3*0.01064/

c
c Input data file formats:
c
c  B*.mrg:
cDay Hr Mn Sec   t-t0 (sec)   BX, BY, BZ, Range
c 666   3  44  12.49        7399.49       77.49979      -66.25665      144.40353 3
c 666   3  54  12.49        7999.49       76.92754      -65.14220      145.40967 3
c
c  2     read(bunit,100,end=99) bday,bhr,bmn,bsec,bt,bx,by,bz,rng
c  100   format(3i4,f7.2,3x,f12.2,3(3x,f12.5),i2)
c
c
c
c
cB-ctp*.mrg:
cDay Hr Mn Sec   t-t0 (sec),BX_ctp,BY_ctp,BZ_ctp
c 666   3  44  12.00        7399.00       -6.01273      -21.33036       28.85257
c 666   3  54  12.00        7999.00       -5.99710      -21.27201       28.77307
c
c        read(bctpunit,110,end=99) bcday,bchr,bcmn,bcsec,bct,bcx,bcy,bcz
c  110   format(3i4,f7.2,3x,f12.2,3(3x,f12.5))
c
c
c
cctp*.mrg:
cDay Hr Mn Sec t-t0 (sec),HTR_INST,HTR_SC  ,HTR_PRP ,ASHNT1  ,DGSHNT  ,SLR_ARRY,STR_CAM ,NM_DPU  ,XG_DPU  ,AIU_1   
c 666   3  44  12.00        7399.00     11      0     56     69    128    160     83     79    177    124
c 666   3  54  12.00        7999.00     17      0     50     70    128    160     83     79    174    124
c
c        read(ctpunit,120,end=99) cday,chr,cmn,csec,ct,inst_htr,sc_htr,prp_htr,ashnt1,dgshnt,slr_arry
c  120   format(3i4,f7.2,3x,f12.2,10i7)
c
c
c  bhsk: added 2/13/00
c
cbhsk*.dat
cDay Hr Mn Sec t-t0 (sec),RATE  ,SPIKES,CAL_ON,OFFSET,RANGE ,MAN_RN,I_DCDC,TPROBE,TELECT,T_DPU 
c 1454   3  17  10.00        3600.00    6.0000    0.0000    0.0000    1.0000    3.0000    3.0000   61.8500   -9.2500   20.3402   15.4634
c 1454   3  17  54.00        3644.00    6.0000    0.0000    0.0000    1.0000    3.0000    3.0000   61.8500   -9.2500   20.3402   15.4634
c 1454   3  18  38.00        3688.00    6.0000    0.0000    0.0000    1.0000    3.0000    3.0000   61.8500   -9.2500   19.7268   15.4634
c 1454   3  19  22.00        3732.00    6.0000    0.0000    0.0000    5.0000    3.0000    3.0000   61.8500   -9.2500   20.3402   15.4634
c 1454   3  20   6.00        3776.00    6.0000    0.0000    0.0000    5.0000    3.0000    3.0000   61.8500   -9.2500   20.3402   15.4634
c

      cos_xy=cosd(thxy)
      sin_xy=sind(thxy)
      cos_xz=cosd(thxz)
      sin_xz=sind(thxz)
      
c
c  Open and read in offsets lookup table:  (added July 7, 2000)
c
c
      open(unit=1,file='offsets.dat',status='old',readonly)
      i=1
7     read(1,*,end=8) xdy0_sc(i),bx0_sc(i),by0_sc(i),bz0_sc(i)
      t0_sc(i)=xdy0_sc(i)*86400.0
      i=i+1
      goto 7
8     nsc0=i-1
      close(unit=1)
      k_sc0=1

c
c  Diagnostic:  (passed 7/7/00)
c
c       do i=1,nsc0
c          type *,xdy0_sc(i),bx0_sc(i),by0_sc(i),bz0_sc(i)
c       enddo

c      type *,'Enter input bYYDDD.mrg file listing:'
c      read(*,50) filelist

      filelist = 'mrg_files.lis'
      
      open(unit=1,file=filelist,status='old',readonly)

c
c  Loop over files in input listing:
c
      
1     read(1,60,end=95) yd
60    format(1x,a5)

      write(bfile,51) yd
      write(outbfile,54) yd
      write(bctpfile,52) yd
      write(ctpfile,53) yd
      write(scfld_file,56) yd
      write(bhskfile,57) yd
      
51    format('b',a5,'.mrg')
52    format('bctp',a5,'.mrg')
53    format('ctp',a5,'.mrg')
54    format('b',a5,'.fix')
56    format('bsc',a5,'.dat')
57    format('bhsk',a5,'.dat')

      bunit=14
      open(unit=bunit,file=bfile,status='old',readonly)
      
      bctpunit=15
      open(unit=bctpunit,file=bctpfile,status='old',readonly)
      
      ctpunit=16
      open(unit=ctpunit,file=ctpfile,status='old',readonly)
c
c  Open bhsk file and read contents of entire file:
c      
      bhsk_unit=18
      open(unit=bhsk_unit,file=bhskfile,status='old',readonly)
      read(bhsk_unit,50) junk

      i=1
21    read(bhsk_unit,*,end=91) dy_hsk,hr_hsk,mn_hsk,s_hsk,dt_hsk,
     &    brate,spikes,cal_on,toffset,rrange,rman_rng,idcdc,tprobe,telect(i),tdpu
      call metd2sec(dy_hsk,hr_hsk,mn_hsk,s_hsk,ttmp)
      tbhsk(i)=ttmp
      i=i+1
      goto 21
      
91    n_hsk=i-1
      k_hsk=1
      close(unit=bhsk_unit)

c
c  Open output SC-field and fix file:
c
      scfld_unit=17
      open(unit=scfld_unit,file=scfld_file,status='unknown')
      
      outbunit=2
      open(unit=outbunit,file=outbfile,status='unknown')
      write(outbunit,55)
55    format('Day Hr Mn Sec   t-t0 (sec)   BX, BY, BZ, Range')
      
c
c  Now loop over entries in the three files:
c

      read(bctpunit,50) junk
      read(bunit,50) junk
      read(ctpunit,50) junk
50    format(a)

      nx=0
      ny=0
      nz=0
      bx_day=0.0
      by_day=0.0
      bz_day=0.0
      today=0

2     read(bunit,*,end=99) bday,bhr,bmn,bsec,bt,bx,by,bz,rng
c100   format(3i4,f7.2,3x,f12.2,3(3x,f12.5),i2)

c
cB-ctp*.mrg:
cDay Hr Mn Sec   t-t0 (sec),BX_ctp,BY_ctp,BZ_ctp
c 666   3  44  12.00        7399.00       -6.01273      -21.33036       28.85257
c 666   3  54  12.00        7999.00       -5.99710      -21.27201       28.77307
c

      read(bctpunit,*,end=99) bcday,bchr,bcmn,bcsec,bct,bcx,bcy,bcz
c110   format(3i4,f7.2,3x,f12.2,3(3x,f12.5))

c
c  Interp telect to tmrg:
c
      call metd2sec(bday,bhr,bmn,bsec,tmrg)
      
      if (tmrg.le.tbhsk(1)) then
         t_interp=telect(1)
     goto 23
      endif
      
      if (tmrg.ge.tbhsk(n_hsk)) then
         t_interp=telect(n_hsk)
     goto 23
      endif
      
22    if (tmrg.gt.tbhsk(k_hsk+1)) then
         k_hsk=k_hsk+1
     if (k_hsk.lt.n_hsk) then
        goto 22
     else
        t_interp=telect(n_hsk)
        goto 23
     endif
      endif
      
24    if (tmrg.le.tbhsk(k_hsk)) then
         k_hsk=k_hsk-1
     if (k_hsk.gt.1) then
        goto 24
     else
        t_interp=telect(1)
        goto 23
     endif
      endif
      
      if ((tmrg.gt.tbhsk(k_hsk)).and.(tmrg.le.tbhsk(k_hsk+1))) then
         t_interp=(   (tmrg-tbhsk(k_hsk))*telect(k_hsk+1) 
     &              + (tbhsk(k_hsk+1)-tmrg)*telect(k_hsk) )
     &              /(tbhsk(k_hsk+1)-tbhsk(k_hsk))
      else
         t_interp=0.5*(tbhsk(1)+tbhsk(n_hsk))
      endif
      
23    continue

c
c  Interp to sc offsets:
c
      if (tmrg.le.t0_sc(1)) then
         boff_SC(1)=bx0_sc(1)
         boff_SC(2)=by0_sc(1)
         boff_SC(3)=bz0_sc(1)
c    type *,'Using first offsets:'
      endif
      
      if (tmrg.ge.t0_sc(nsc0)) then
         boff_SC(1)=bx0_sc(nsc0)
         boff_SC(2)=by0_sc(nsc0)
         boff_SC(3)=bz0_sc(nsc0)
c    type *,'Using last offsets:',boff_SC(1),boff_SC(2),boff_SC(3)
      endif
      
      if ((tmrg.gt.t0_sc(1)).and.(tmrg.lt.t0_sc(nsc0))) then
      
31       if (tmrg.gt.t0_sc(k_sc0+1)) then
        k_sc0=k_sc0+1
        type *,'k_sc0 changed to: ',k_sc0
        goto 31
     endif
     if (tmrg.le.t0_sc(k_sc0)) then
        k_sc0=k_sc0-1
        type *,'k_sc0 changed to: ',k_sc0
        goto 31
     endif
     
     wt0=(t0_sc(k_sc0+1)-tmrg)/(t0_sc(k_sc0+1)-t0_sc(k_sc0))
     wt1=(tmrg-t0_sc(k_sc0))/(t0_sc(k_sc0+1)-t0_sc(k_sc0))
     boff_SC(1) = wt0*bx0_sc(k_sc0) + wt1*bx0_sc(k_sc0+1)
     boff_SC(2) = wt0*by0_sc(k_sc0) + wt1*by0_sc(k_sc0+1)
     boff_SC(3) = wt0*bz0_sc(k_sc0) + wt1*bz0_sc(k_sc0+1)
c    type *,'Using interpolated offsts:',k_sc0,boff_SC(1),boff_SC(2),boff_SC(3)
     
      endif

c
cctp*.mrg:
cDay Hr Mn Sec t-t0 (sec),HTR_INST,HTR_SC  ,HTR_PRP ,ASHNT1  ,DGSHNT  ,SLR_ARRY,STR_CAM ,NM_DPU  ,XG_DPU  ,AIU_1   ,TubeTemp,ProbTemp   
c 666   3  44  12.00        7399.00     11      0     56     69    128    160     83     79    177    124 -5.14  -0.45
c 666   3  54  12.00        7999.00     17      0     50     70    128    160     83     79    174    124 -5.14  -0.45
c

      read(ctpunit,*,end=99) cday,chr,cmn,csec,ct,inst_htr,sc_htr,prp_htr,ashnt1,dgshnt,slr_arry,
     &                       str_cam,nm_dpu,xg_dpu,aiu_1,tube_temp,probe_temp
c120   format(3i4,f7.2,3x,f12.2,10i7)

      dig_rat=float(dgshnt)/float(slr_arry)
      call get_dg_index(dig_rat,id,dg_lowr,dg_uppr,20)
      
c
c  B-ctp, analog shunt, digital shunt step, DC offset and solar array current
c  corrections:
c


c -----------------Version 1 Original------------------------
c      
c      bx=bx
c     &     -bcx                          ! Terminal board
c     &     -bx_ashnt(id)*float(ashnt1)   ! Analog shunt
c     &     -bx_dshunt(id)                ! Digital shunt rel. to #18
c     &     -arry_bx*float(slr_arry)      ! Solar array total current
c     
c      by=by
c     &     -bcy                          ! Terminal board
c     &     -by_ashnt(id)*float(ashnt1)   ! Analog shunt
c     &     -by_dshunt(id)                ! Digital shunt rel. to #18
c     &     -arry_by*float(slr_arry)      ! Solar array total current
c     
c      bz=bz
c     &     -bcz                          ! Terminal board
c     &     -bz_ashnt(id)*float(ashnt1)   ! Analog shunt
c     &     -bz_dshunt(id)                ! Digital shunt rel. to #18
c     &     -arry_bz*float(slr_arry)      ! Solar array total current
c-------------------------------------------------------------


c -----------------Version 2  12/14/98-------------------------
c      
c      bx=bx
c     &     -bcx                          ! Terminal board
c     &     -xash(ashnt1)            ! Analog shunt
c     &     -bx_dshunt(id)                ! Digital shunt rel. to #0
c     &     -arry_bx*float(slr_arry)      ! Solar array total current
c      
c      by=by
c     &     -bcy                          ! Terminal board
c     &     -yash(ashnt1)            ! Analog shunt
c     &     -by_dshunt(id)                ! Digital shunt rel. to #0
c     &     -arry_by*float(slr_arry)      ! Solar array total current
c      
c      bz=bz
c     &     -bcz                          ! Terminal board
c     &     -zash(ashnt1)            ! Analog shunt
c     &     -bz_dshunt(id)                ! Digital shunt rel. to #0
c     &     -arry_bz*float(slr_arry)      ! Solar array total current
c-------------------------------------------------------------


c -----------------Version 3b  12/28/98-------------------------
c
c      bash(1)=xash(ashnt1)
c      barry(1)=arry_bx(id)*float(slr_arry)
c      
c      bxSC=bcx                           ! Terminal board
c     &     +bash(1)                     ! Analog shunt
c     &     +bx_dshunt(id)                ! Digital shunt rel. to #0
c     &     +barry(1)                     ! Solar array total current
c     
c      bx=bx-bxSC
c
c      bash(2)=yash(ashnt1)
c      barry(2)=arry_by(id)*float(slr_arry)
c      
c      bySC=bcy                           ! Terminal board
c     &     +bash(2)                     ! Analog shunt
c     &     +by_dshunt(id)                ! Digital shunt rel. to #0
c     &     +barry(2)                     ! Solar array total current
c     
c      by=by-bySC
c
c      bash(3)=zash(ashnt1)
c      barry(3)=arry_bz(id)*float(slr_arry)
c      
c      bzSC=bcz                           ! Terminal board
c     &     +bash(3)                     ! Analog shunt
c     &     +bz_dshunt(id)                ! Digital shunt rel. to #0
c     &     +barry(3)                     ! Solar array total current
c     
c      bz=bz-bzSC
c      
c-------------------------------------------------------------

c -----------------Version 4  1/4/99-------------------------
c
c      bash(1)=xash(ashnt1)
c      barry(1)=arry_bx(id)*float(slr_arry)
c      bdshunt(1)=bx_dshunt(id)*float(slr_arry)
c      btelec(1)=2.17 - 0.0889*t_interp
c      
c      bxSC=bcx                           ! Terminal board
c     &     +bash(1)                     ! Analog shunt
c     &     +bdshunt(1)                   ! Digital shunt rel. to #0
c     &     +barry(1)                     ! Solar array total current
c     &     +btelec(1)                    ! Temperature dependent offset.
c     
c      bx=bx-bxSC
c
c      bash(2)=yash(ashnt1)
c      barry(2)=arry_by(id)*float(slr_arry)
c      bdshunt(2)=by_dshunt(id)*float(slr_arry)
c      btelec(2)=-0.10 - 0.1347*t_interp
c      
c      bySC=bcy                           ! Terminal board
c     &     +bash(2)                     ! Analog shunt
c     &     +bdshunt(2)                   ! Digital shunt rel. to #0
c     &     +barry(2)                     ! Solar array total current
c     &     +btelec(2)                    ! Temperature dependent offset.
c     
c      by=by-bySC
c
c      bash(3)=zash(ashnt1)
c      barry(3)=arry_bz(id)*float(slr_arry)
c      bdshunt(3)=bz_dshunt(id)*float(slr_arry)
c      btelec(3)=3.12 - 0.1352*t_interp
c      
c      bzSC=bcz                           ! Terminal board
c     &     +bash(3)                     ! Analog shunt
c     &     +bdshunt(3)                   ! Digital shunt rel. to #0
c     &     +barry(3)                     ! Solar array total current
c     &     +btelec(3)                    ! Temperature dependent offset.
c     
c      bz=bz-bzSC
c      
c-------------------------------------------------------------

c -----------------Version 5  7/7/00-------------------------

      bash(1)=xash(ashnt1)
      barry(1)=arry_bx(id)*float(slr_arry)
      bdshunt(1)=bx_dshunt(id)*float(slr_arry)
      btelec(1)=2.17 - 0.0889*t_interp
      
      bxSC=bcx                           ! Terminal board
     &     +bash(1)                      ! Analog shunt
     &     +bdshunt(1)                   ! Digital shunt rel. to #0
     &     +barry(1)                     ! Solar array total current
     &     +btelec(1)                    ! Temperature dependent offset.
     &     +boff_SC(1)                   ! 27-day time scale variable SC offset.
     
      bx=bx-bxSC

      bash(2)=yash(ashnt1)
      barry(2)=arry_by(id)*float(slr_arry)
      bdshunt(2)=by_dshunt(id)*float(slr_arry)
      btelec(2)=-0.10 - 0.1347*t_interp
      
      bySC=bcy                           ! Terminal board
     &     +bash(2)                      ! Analog shunt
     &     +bdshunt(2)                   ! Digital shunt rel. to #0
     &     +barry(2)                     ! Solar array total current
     &     +btelec(2)                    ! Temperature dependent offset.
     &     +boff_SC(2)                   ! 27-day time scale variable SC offset.
     
      by=by-bySC

      bash(3)=zash(ashnt1)
      barry(3)=arry_bz(id)*float(slr_arry)
      bdshunt(3)=bz_dshunt(id)*float(slr_arry)
      btelec(3)=3.12 - 0.1352*t_interp
      
      bzSC=bcz                           ! Terminal board
     &     +bash(3)                      ! Analog shunt
     &     +bdshunt(3)                   ! Digital shunt rel. to #0
     &     +barry(3)                     ! Solar array total current
     &     +btelec(3)                    ! Temperature dependent offset.
     &     +boff_SC(3)                   ! 27-day time scale variable SC offset.
     
      bz=bz-bzSC
      
c
c  Tank valve open correction.  Applies after 157414531s met.
c
c    bxS=53.22
c    byS=-8.96
c    bzS=64.69
c    if (tmet Ñ 157414531) then
c       bx= bx-bxS
c       by=by-byS
c       bz=bz-bzS
c    endif
c
      if (tmrg.ge.157414531.0) then
         bx=bx-53.22
     by=by+8.96
     bz=bz-64.69
      endif
      
c-------------------------------------------------------------

c
c  I have no idea why this was restricted with respect to the met day.
c  575 was the day that Tamitha Mulligan and I analyzed with respect to 
c  offsets.  This was the checkout after the upload.
c  The software upload was day 572.  Day 580 was the day the automomy
c  rule for CTP-MAG synchronization was implemented.  The offset correction
c  should always be made.
c
c  Besides this is done in getdata.pro.  Bizzare.
c
c      if (bday.lt.578) then
c         bx=bx-ox(rng)
c         by=by-oy(rng)
c         bz=bz-oz(rng)
c      endif
c
c  I suspect there was a sign error in getdata.pro.  To check for this
c  I temporarily add 2*ox(rng) back into the data.  This is for diagnostic
c  purposes only!  2/13/97
c

c      type *,'rng, ox,oy,oz:',rng,ox(rng+1),oy(rng+1),oz(rng+1)

c      bx=bx+2.0*ox(rng+1)
c      by=by+2.0*oy(rng+1)
c      bz=bz+2.0*oz(rng+1)
      
c      if (bday.lt.300) then
c         bx=bx-bxoffset(1)
c         by=by-byoffset(1)
c         bz=bz-bzoffset(1)
c      endif
      
c      if (bday.ge.300) then
c         bx=bx-bxoffset(2)
c         by=by-byoffset(2)
c         bz=bz-bzoffset(2)
c      endif

c   Version 3b 12/28/98

      if (bday.lt.300) then
         boffset(1)=bxoffset(1)
         boffset(2)=byoffset(1)
         boffset(3)=bzoffset(1)
      endif
      
      if (bday.ge.300) then
         boffset(1)=bxoffset(2)
         boffset(2)=byoffset(2)
         boffset(3)=bzoffset(2)
      endif
      
      bx=bx-boffset(1)
      by=by-boffset(2)
      bz=bz-boffset(3)
      
c
c  Apply sensor orientation correction: 2/20/98
c  Determined from Earth flyby.
c

      bxtmp=cos_xy*bx+sin_xy*by    ! 1.76 deg rotation in XY-plane
      bytmp=cos_xy*by-sin_xy*bx
      
      bx=bxtmp
      by=bytmp
      
      bxtmp=cos_xz*bx+sin_xz*bz    ! 0.22 deg rotation in X'Z-plane
      bztmp=cos_xz*bz-sin_xz*bx
      bx=bxtmp
      bz=bztmp
     
      
      write(outbunit,101) bday,bhr,bmn,bsec,bt,bx,by,bz,rng
101   format(i5,2i4,f7.2,3x,f12.2,3(3x,f14.5),i2)

      write(scfld_unit,102) bday,bhr,bmn,bsec,bt,bxSC,bySC,bzSC,id,
     &  bcx,bcy,bcz,(bash(i),i=1,3),
     &  (barry(i),i=1,3),(bdshunt(i),i=1,3),(btelec(i),i=1,3)
102   format(i5,2i4,f7.2,3x,f12.2,3(3x,f14.5),1x,i2,15f8.2)
      
      goto 2


99    close(unit=bunit)
      close(unit=bctpunit)
      close(unit=ctpunit)
      
      write(*,15) yd
15    format(1x,'Day ',a5,' done.')


      goto 1
      

95    close(unit=1) 
      
      stop
      end
c
c-------------------------------------------------------------------------
c
      subroutine get_dg_index(dig_rat,k,dg_lowr,dg_uppr,nmax)
      real dig_rat,dg_uppr(20),dg_lowr(20)
      integer k,i,nmax
      
      i=0
      k=0
      do i=1,nmax
         if ((dig_rat.ge.dg_lowr(i)).and.(dig_rat.le.dg_uppr(i))) k=i
      enddo
      
      return
      end
c
c-------------------------------------------------------------------------
c
      real function xash(va)
      real xtmp,wd,wt0,wt1,xpnts(7)
      integer va,i,ab(7)
    
      data ab,xpnts/0,17,34,51,68,85,102,
     &  -3.18863, -2.14799, -1.55278, -0.1651,
     &   1.3744,   1.8969,   2.1785/
    
      i=2
1     if (va.gt.ab(i)) then
         i=i+1
     if (i.gt.7) then
            xtmp=xpnts(7)
        goto 2
     endif
     goto 1
      endif
      
    wd=float(ab(i)-ab(i-1))
    wt0=float(ab(i)-va)
    wt1=float(va-ab(i-1))
    xtmp=(wt0*xpnts(i-1)+wt1*xpnts(i))/wd
c      type *,va,ab(i-1),ab(i),xpnts(i-1),xpnts(i),xtmp
    
2      xash=xtmp
      end
c
c-------------------------------------------------------------------------
c
      real function yash(va)
      real ytmp,wd,wt0,wt1,ypnts(7)
      integer va,i,ab(7)
    
      data ab,ypnts/0,17,34,51,68,85,102,
     &   -2.7746, -1.45173, -0.9096, -0.06753,
     &    0.64605,  2.5994,  3.4384/
    
      i=2
1     if (va.gt.ab(i)) then
         i=i+1
     if (i.gt.7) then
            ytmp=ypnts(7)
        goto 2
     endif
     goto 1
      endif
      
    wd=float(ab(i)-ab(i-1))
    wt0=float(ab(i)-va)
    wt1=float(va-ab(i-1))
    ytmp=(wt0*ypnts(i-1)+wt1*ypnts(i))/wd
c      type *,va,ab(i-1),ab(i),ypnts(i-1),ypnts(i),ytmp
    
2     yash=ytmp
      end
c
c-------------------------------------------------------------------------
c
      real function zash(va)
      real ztmp,wd,wt0,wt1,zpnts(7)
      integer va,i,ab(7)
    
      data ab,zpnts/0,17,34,51,68,85,102,
     &   -0.1554, -0.2329,  0.5596, -0.3298,
     &   -0.1260,  0.1066, -0.7398/
    
      i=2
1     if (va.gt.ab(i)) then
         i=i+1
     if (i.gt.7) then
            ztmp=zpnts(7)
        goto 2
     endif
     goto 1
      endif
      
      
    wd=float(ab(i)-ab(i-1))
    wt0=float(ab(i)-va)
    wt1=float(va-ab(i-1))
    ztmp=(wt0*zpnts(i-1)+wt1*zpnts(i))/wd
    
2     zash=ztmp
      end
c
c-------------------------------------------------------------------------
c
      subroutine metd2sec(d,h,m,s,tsec)
      integer d,h,m
      real s
      real*8 tsec
      
      tsec=86400.0*dfloat(d) + 3600.0*dfloat(h) + 60.0*dfloat(m) + s
      return
      end
c
c-------------------------------------------------------------------------
c


--------------------------
 530.000   -1.118   -0.071    0.197
 557.957   -1.118   -0.071    0.197
 587.972   -0.518    0.077   -0.210
 616.040   -1.853    0.259    0.818
 643.417   -1.980   -0.069    1.207
 672.808   -1.653    0.001   -0.960
 702.683   -3.804    0.667   -0.860
 750.000   -3.804    0.667   -0.860
 751.000    0.215   -0.110   -2.593
 894.436    0.215   -0.110   -2.593
 921.684    0.458    0.336   -1.247
 951.705    0.618    0.678   -1.393
 981.634    1.510   -0.079    0.546
1011.374   -0.260   -0.726    2.348
1036.070    0.643   -1.675    1.074
1070.000    0.643   -1.675    1.074
1071.000    0.997    2.491    0.992
1095.009    0.997    2.491    0.992
1105.000    0.997    2.491    0.992
1106.000    1.189    0.163    3.568
1144.313    1.189    0.163    3.568
1172.639   -0.655   -0.233    1.680
1200.307   -0.172   -0.293    1.200
1229.570    0.366   -1.144   -0.858
1258.983    0.896   -0.950    1.148
1287.725    0.597   -0.540   -1.093
1315.285    0.387   -0.719    0.509
1342.416   -0.516   -0.981    0.729
1369.470    0.985   -1.765    0.371
1396.845    0.414   -0.992    0.539
1424.394    0.617   -1.676    0.347
1453.839   -0.034    0.312    0.888
1482.086   -0.223   -0.958    2.241
1509.480    0.883   -1.861    2.669
1536.955    1.119   -1.270    2.733
1564.121    0.689   -1.427    3.328
1591.244    1.326   -2.163    2.524
1618.345    0.013   -3.027    1.780
1645.418    1.157   -2.446    2.295
1673.249    0.007   -4.572    3.224
1701.569    1.352   -1.345    1.943
1729.625    1.167   -1.080    1.288
1757.285    0.661   -1.855    1.995
1791.715    0.310   -0.275    1.710
1813.867    0.182   -0.353    2.950
2000.000    0.182   -0.353    2.950
$Name: $
$Id: $
1 /project/near/SDC_archive/SpiceOps/ops/spk/near_1996151_2000051_merged_v01.bsp
1 /project/near/SDC_archive/SpiceOps/ops/spk/near_1999218_2000101_merged_v02.bsp
1 /project/near/SDC_archive/SpiceOps/ops/spk/near_2000040_2000125_reconstr_v02.bsp
1 /project/near/SDC_archive/SpiceOps/ops/spk/near_2000166_2000291_merged_v01.bsp
1 /project/near/SDC_archive/SpiceOps/ops/spk/near_2000121_2000215_reconstr_v01.bsp
1 /project/near/SDC_archive/SpiceOps/ops/spk/near_2000222_2001028_merged_v03.bsp
1 /project/near/SDC_archive/SpiceOps/ops/spk/near_2000300_2001043_merged_v03.bsp
1 /project/near/SDC_archive/SpiceOps/ops/spk/near_2001020_2001060_reconstr_v01.bsp
0 /project/near/SDC_archive/SpiceOps/ops/current/near.tscl
0 /project/near/SDC_archive/SpiceOps/ops/current/near.tls
0 /project/near/SDC_archive/SpiceOps/ops/ik/grs12.ti
0 /project/near/SDC_archive/SpiceOps/ops/ik/msi15.ti
0 /project/near/SDC_archive/SpiceOps/ops/ik/nis14.ti
0 /project/near/SDC_archive/SpiceOps/ops/ik/nlr12.ti
0 /project/near/SDC_archive/SpiceOps/ops/ik/xrs12.ti
0 /project/near/SDC_archive/SpiceOps/ops/pck/pck00007.tpc
2 /project/near/SDC_archive/SpiceOps/ops/pck/earth_970804_200101_pred.bpc
2 /project/near/SDC_archive/SpiceOps/ops/pck/erosatt_1998329_2001157_v01.bpc
0 /project/near/SDC_archive/SpiceOps/ops/pck/eros_1999008.tpc
108