Program Flat(Input, Output) ;
{ Reads all specified images, a strip at a time }
{ Forms an image equal to the median of scaled deviations of }
{ pixels from a precomputed median }
Uses WinAPI, X_U, FITSX_U, SORT_U ;
Const
  DarkName   = 'C:\BP\TASS\CD5\H4RDARK.FTS' ; { Input Dark file }
  FlatName   = 'C:\BP\TASS\CD5\H4RFLAT.FTS' ; { Output Flat file }
  MedianName = 'C:\BP\TASS\CD5\H4RMED.TXT' ;  { Input file of medians }
{ Analysis area }
  Min2 : Integer   =    2 ;
  Max2 : Integer   = 2033 ;
  Min1 : Integer   =    6 ;
  Max1 : Integer   = 2037 ;
  NFil = 35 ;
  Name : Array[0..NFil-1] of String =
   ('D:\H4R1438.878',
    'D:\H4R1438.882',
    'D:\H4R1438.886',
    'D:\H4R1438.890',
    'D:\H4R1438.895',
    'D:\H4R1438.899',
    'D:\H4R1438.903',
    'D:\H4R1439.870',
    'D:\H4R1439.875',
    'D:\H4R1439.880',
    'D:\H4R1439.885',
    'D:\H4R1439.891',
    'D:\H4R1439.896',
    'D:\H4R1439.901',
    'D:\H4R1440.871',
    'D:\H4R1440.875',
    'D:\H4R1440.879',
    'D:\H4R1440.884',
    'D:\H4R1440.888',
    'D:\H4R1440.892',
    'D:\H4R1440.897',
    'D:\H4R1442.863',
    'D:\H4R1442.867',
    'D:\H4R1442.871',
    'D:\H4R1442.876',
    'D:\H4R1442.880',
    'D:\H4R1442.885',
    'D:\H4R1442.889',
    'D:\H4R1446.853',
    'D:\H4R1446.858',
    'D:\H4R1446.862',
    'D:\H4R1446.866',
    'D:\H4R1446.871',
    'D:\H4R1446.875',
    'D:\H4R1446.880') ;
Var
  F : Array[0..NFil-1] of FITSInfoT ;
  Drk : FITSInfoT ;
  Flt : FITSInfoT ;
  FNo : Integer ;
  Med : Array[0..NFil-1] Of Longint ; { Longint to thwart BP7.0 compiler }
  MedDrk : Longint ;
  I2, I1, N2_M, N1_M : Longint ;
  L : Array[0..NFil-1] of Single ; { Started out as L for Longint. Um! }
  Index : Array[0..NFil-1] of Integer ;
  MedFile : Text ;
  W : Longint ;{ Width of strip }
  S : String ;
  Val  : Longint ;
  IVal : Integer ;
  T : Single ;
Begin
  N2_M := Max2 - Min2 + 1 ;
  N1_M := Max1 - Min1 + 1 ;
  FileMode := 0 ; { Read only }
  Assign(MedFile, MedianName) ;
  Reset(MedFile) ;
  Readln(MedFile, S) ; { Throw away header without looking at it }
{ Nobody would be silly enough to supply the wrong file!         }
  Readln(MedFile, MedDrk) ;
  For FNo := 0 To NFil-1 Do Readln(MedFile, Med[FNo]) ;
  Close(MedFile) ;
  SetHead(Flt, N1_M, N2_M) ;
  W := 64 ; { Width of strip to be read each time }
  Repeat
    If (F[1].Off2 + W > Max2) Then W := Max2 - F[1].Off2 + 1 ;
{ This is a bit Byzantine: ReadStrip and WriteStrip somehow know }
{ when they are being called for the first time on a particular  }
{ file and act accordingly. I must really get around to fixing   }
{ them so that this is made clear!                               }
    ReadStrip(Drk, DarkName, 0, N1_M, 0, W) ;
    For FNo := 0 to NFil-1 Do Begin
      ReadStrip(F[FNo], Name[FNo], Min1, N1_M, Min2, W) ;
    End ;
    Writeln('Processing') ;
    For I2 := 0 to W-1 Do For I1 := 0 to N1_M-1 Do Begin
      For FNo := 0 to NFil-1 Do Begin
{ The following contortions get round BP7.0 compiler bugs }
        T := IPtr(F[FNo].Image, I2*N1_M+I1)^ ; { IPtr: 32-bit adressing }
{ equivalent to F.Image[I2, I1]                                         }
        L[FNo] := (T - IPtr(Drk.Image, I2*N1_M+I1)^)/(Med[FNo] - MedDrk) ;
      End ;
      SSort(L, 0, NFil-1, Index) ; { SSort = Single prec FP QuickSort }
      T := L[Index[(NFil-1) Div 2]] ;
{ Compute the flat field scale factor. Note this is the reciprocal of }
{ the MKIII program convention.                                       }
      If (T > 0.1) Then Val := Round(32768/T) Else Val := 327680 ;
      If (Val > 65534) Then Val := 65534 ; { Fireproof }
      If (Val < 16384) Then Val := 16384 ;
      IVal := Val - 32768 ;
      IPtr(Flt.Image, I2*N1_M+I1)^ := IVal ;
    End ;
    Writeln('Writing from ', Flt.Off2, '  ', W) ;
    WriteStrip(Flt, FlatName, 99, N1_M, 0, W) ; { 99 is a flag }
{ See notes above re: Byzantine }
    Writeln('Written to   ', Flt.Off2) ;
  Until (Flt.Off2 >= N2_M) ;
End. { Flat }