	subroutine makefile(ch)
c
C	THIS IS MAKEFILE20
C
C	A PROGRAM TO SEARCH FOR LANGMUIR WAVES, MOSTLY UPSTREAM
C		 AND TO SAVE AMPLITUDE AND AMPLITUDE OF HARMONIC
C
C	ORIGINALLY, THE PROGRAM DID NOT COMPUTE ANGLE TO B.  ON 13 SEP 2007 
C	THIS WAS ADDED AND A NEW SET OF RESULTS WRITTEN TO FOR080.DAT
C	note, 9 June 2008  I find only two criteria, bandwidth lt .02 freq
C		and freq between .75 fp and 1.25 fp
C
	COMMON /TDSF_BACKGROUND/ TDSFBCK(1025,2,4),NBCK(2,4)
C
	integer*4 ch,ok,okt,OK1,OK2,SCETI4(2),N1DATA(2048),N2DATA(2048)
	INTEGER*4 W_CHANNEL_CLOSE,W_EVENT,RET_SIZE,W_MESSAGES_OFF
	INTEGER*4 W_ITEM_I4,W_ITEM_R4,W_ITEM_R8
	INTEGER*4 TDS_CHANNEL,ISRC,SUNCLOCK,ERTDAY,ERT(2)
	REAL*8 SCET8,XGSE,YGSE,ZGSE,RGSE
	REAL DBSPEC1(1025),V1DATA(2048)
	REAL DBSPEC2(1025),V2DATA(2048)
	REAL ZCROSS(1024),Z1INT(1024),Z2INT(1024)
	REAL FREQ(1025),FCOUNT(2),FAVR(2),FSTD(2),F3MOM(2),FBW(2),SKEW(2)
	REAL F2COUNT(2),F2AVR(2),F2STD(2),F2BW(2)
	REAL FPK(5),BWPK(5),AMPK(5)
	REAL VMATRIX(2,2),EVAL(2)
	character*32 ITEM
	character*4 event
	DATA TWOPI /6.2831853/
	DATA YXRATIO /11.8/
C	ABOVE BASED ON LEFFX = 41.1, LEFFY = 3.79
	DATA PXTOT,PYTOT /0.,0./
	DATA RE /6.378E3/
	DATA ANGTOB,TDIST /0.,0./
	DATA NFLLIMP,NFUPLIMP /40,1023/
	data event /'TDSF'/
C	data event /'FILL'/
C
	YXSQ = YXRATIO**2
	IFOUND = 0
C
C
C	CALCULATE FREQUENCIES FOR SPECTRUM
C
	SPS = 120000.
	DO N = 1,1025
	  FREQ(N) = (N-1)*SPS/2048.    ! FREQ(N) GOES WITH DBSPECT(N)
C	  DBSPECT(2) IS LOWEST NON-ZERO FREQUENCY, DBSPECT(1) IS DC
	ENDDO
C
C	OPEN(UNIT=90,FILE='MAKEFILE20.RESULTS',STATUS='OLD',ACCESS='APPEND')
 100	ok = w_event(ch,event)
C
C	CHECK FOR END OF RECORD
C
	if (ok.ne.1) then
		if(ok.eq.82) then
		   ok = w_channel_close(ch)
		   CLOSE(UNIT=90)
		   return
	        endif
		write(6,*) 'cannot open ',event, ', ok=', ok
	endif
C
C	END OF END OF RECORD CHECK
C
	OKT = W_MESSAGES_OFF(ch)
C
C
	ITEM = 'CHANNEL'
	ok = W_ITEM_I4(ch, item, TDS_CHANNEL, 1, ret_size)
C
	IPROS = 4
	IF(TDS_CHANNEL.EQ.1) THEN
C	  ITEM = 'DATA'
C	  ok1 = W_ITEM_I4(ch, item, N1DATA, 2048, ret_size)
	  CALL TDS_PHYS(CH,IPROS,N1DATA,V1DATA,DBSPEC1)
	  ITEM = 'EVENT_NUMBER'
	  ok = W_ITEM_I4(ch, item, NO_EVT1, 1, ret_size)
		MAXCH1 = 0
		DO N = 1,2048
		  MAXCH1 = MAX0(MAXCH1,IABS(N1DATA(N)-128)) 
		ENDDO
	  GO TO 100
	ELSEIF(TDS_CHANNEL.EQ.2) THEN
C	  ITEM = 'DATA'
C	  ok2 = W_ITEM_I4(ch, item, N2DATA, 2048, ret_size)
	  CALL TDS_PHYS(CH,IPROS,N2DATA,V2DATA,DBSPEC2)
	  ITEM = 'EVENT_NUMBER'
	  ok = W_ITEM_I4(ch, item, NO_EVT, 1, ret_size)
	  IF(NO_EVT.NE.NO_EVT1) GO TO 100
		MAXCH2 = 0
		DO N = 1,2048
		  MAXCH2 = MAX0(MAXCH2,IABS(N2DATA(N)-128)) 
		ENDDO
	ELSE
	  GO TO 100
	ENDIF

C
CDIAG	PRINT*,'******START EVENT',NO_EVT1,NO_EVT,' *****'
	PRINT*,'******START EVENT',NO_EVT1,NO_EVT,' *****'
C
	IF(TDS_CHANNEL.EQ.1) GO TO 100		! GO GET CHANNEL 2
	IF(TDS_CHANNEL.NE.2) GO TO 100
C
C	CALCULATE VMATRIX, FOR ANGLE TO B
C
        VMATRIX(1,1) = 0.
        VMATRIX(1,2) = 0.
        VMATRIX(2,2) = 0.
	XEFFLEN = 41.1
	YEFFLEN = 3.79
        DO N = 1,2048
          VMATRIX(1,1) = VMATRIX(1,1) + (V1DATA(N)/XEFFLEN)**2
          VMATRIX(1,2) = VMATRIX(1,2) + V1DATA(N)*V2DATA(N)/XEFFLEN/
     1	YEFFLEN
          VMATRIX(2,2) = VMATRIX(2,2) + (V2DATA(N)/YEFFLEN)**2
        ENDDO 
C
	XRMS = 1000.*SQRT(VMATRIX(1,1)/2048.)    ! in mV/m
	YRMS = 1000.*SQRT(VMATRIX(2,2)/2048.)
C
	ITEM = 'WIND_ORBIT_X(GSE)_R8'
	ok = W_ITEM_R8(ch, item, XGSE, 1, ret_size)
	XRE = XGSE/RE
	ITEM = 'WIND_ORBIT_Y(GSE)_R8'
	ok = W_ITEM_R8(ch, item, YGSE, 1, ret_size)
	YRE = YGSE/RE
	ITEM = 'WIND_ORBIT_Z(GSE)_R8'
	ok = W_ITEM_R8(ch, item, ZGSE, 1, ret_size)
	ZRE = ZGSE/RE
	RRE = SQRT(XRE**2+YRE**2+ZRE**2)
C
C	IF(RGSE.LT.6.4D05) THEN			! LIMIT 100 RE
C	    ADVANCE BY ONE HOUR AND GET ANOTHER EVENT
C	    ITEM = 'EVENT_SCET_R8'
C	    OK = W_ITEM_R8(ch, item, SCET8, 1, ret_size)
	    ITEM = 'EVENT_SCET'
	    ok = W_ITEM_I4(ch, item, SCETI4, 2, ret_size)
C	    PRINT*,RGSE/6.38E3,'INSIDE 100 RE',SCETI4
C	    SCET8 = SCET8 + 1.D00/24.D00
C	    CALL W_CHANNEL_POSITION(CH,SCET8)
C	    GO TO 100
C	  ENDIF
C
	ITEM = 'SOURCE'
	ok = W_ITEM_I4(ch, item, ISRC, 1, ret_size)
	ITEM = 'EVENT_NUMBER'
	ok = W_ITEM_I4(ch, item, NO_EVT, 1, ret_size)
	ITEM = 'EVENT_SCET'
	ok = W_ITEM_I4(ch, item, SCETI4, 2, ret_size)
	ITEM = 'FAST_RX_SPEED'
	ok = W_ITEM_I4(ch, item, ISPS, 1, ret_size)
	IF(ISPS.NE.0) GO TO 100
	ITEM = 'WIND_MFI_BX(GSE)_R4'
	ok = W_ITEM_R4(ch, item, BX, 1, ret_size)
	ITEM = 'WIND_MFI_BY(GSE)_R4'
	ok = W_ITEM_R4(ch, item, BY, 1, ret_size)
	ITEM = 'WIND_MFI_BZ(GSE)_R4'
	ok = W_ITEM_R4(ch, item, BZ, 1, ret_size)
	ITEM = 'WIND_MFI_BPHI(GSE)_R4'
	OK = W_ITEM_R4(CH, ITEM, AZMAG,1,RET_SIZE)
	ITEM = 'WIND_MFI_BTHETA(GSE)_R4'
	OK = W_ITEM_R4(CH, ITEM, THETAB,1,RET_SIZE)
	ITEM = 'SUN_ANGLE'
	ok = W_ITEM_I4(ch, item, SUNCLOCK, 1, ret_size)
C	ITEM = 'WIND_3DP_E_TEMP_R4'
C	ok = W_ITEM_R4(ch, item, TEMPE, 1, ret_size)
C	ITEM = 'WIND_3DP_ION_TEMP_R4'
C	ok = W_ITEM_R4(ch, item, TEMPI, 1, ret_size)
	ITEM = 'WIND_3DP_ION_DENSITY_R4'
	ok = W_ITEM_R4(ch, item, DENS, 1, ret_size)
C	ITEM = 'WIND_3DP_ION_VX(GSE)_R4'
C	ok = W_ITEM_R4(ch, item, VX, 1, ret_size)
C	ITEM = 'WIND_3DP_ION_VY(GSE)_R4'
C	ok = W_ITEM_R4(ch, item, VY, 1, ret_size)
C	ITEM = 'WIND_MFI_BMAG_R4'
C	ok = W_ITEM_R4(ch, item, BMAG, 1, ret_size)
	ITEM = 'WIND_SPIN_RATE_R4'
	ok = W_ITEM_R4(ch, item, SPINRATE, 1, ret_size)
	if(ok.ne.1) then
	  spinrate = 2.
	endif
	ITEM = 'EVENT_TM_SCET_I4'
	ok = W_ITEM_I4(ch, item, ERT, 2, ret_size)
	ERTDAY = MOD(ERT(1),100)
C
	EANGLE =  -360.*(SUNCLOCK-14.)/4096. - 45. ! ANGLE SUN TO +EX AT END
	IF(EANGLE.LT.-180.) EANGLE = EANGLE + 360.
	IF(EANGLE.GT.180.)  EANGLE = EANGLE - 360.
	DANG = SPINRATE*360./SPS/TWOPI
	THANT = EANGLE + 3072.*DANG  ! ANGLE SUN TO +EX AT START 16nov99 
	THANTM = THANT + 1024.*DANG  ! ANGLE SUN TO +EX AT MIDDLE 16nov99 
C
C	EVALUATE MINIMUM VARIANCE STUFF USING FACT THAT TRACE
C		AND DETERMINANT ARE UNCHANGED BY ROTATIONS.
C		THIS IS ENOUGH FOR A 2 X 2 MATRIX
C
	DET = VMATRIX(1,1)*VMATRIX(2,2) - VMATRIX(1,2)**2
	TRACE = VMATRIX(1,1) + VMATRIX(2,2)
	DISC = SQRT(TRACE**2 - 4.*DET)
	EVAL(1) = .5*(TRACE + DISC)
	EVAL(2) = .5*(TRACE - DISC)
	AXRATIO = SQRT(EVAL(2)/EVAL(1))
C	AX MAJOR AND AX MINOR ARE RATIOS OF Y COMPONENT TO X COMPONENT OF
C		SLOPES OF MAJOR AND MINOR EIGENVECTORS
C		AND GIVE DIRECTION OF E.
c
	AXMAJOR = (EVAL(1)-VMATRIX(1,1))/VMATRIX(1,2)
	AXMINOR = (EVAL(2)-VMATRIX(1,1))/VMATRIX(1,2)
C	SLOPE = YXRATIO*AXMAJOR
	SLOPE = AXMAJOR
	print*,'vmatrix',vmatrix(1,1),vmatrix(1,2),vmatrix(2,2)
	GSEANG = THANTM - ATAND(SLOPE)
	ANGTOB = GSEANG - AZMAG
c	print*,'eval',eval
c	print*,'thant',thant
c	print*,'azmag,bx,by',azmag,bx,by
c	print*,'axmajor,axminor,slope',axmajor,axminor,slope
c	print*,'angtob',angtob
	IF(ANGTOB.GT.90.) ANGTOB = ANGTOB - 180.
	IF(ANGTOB.GT.90.) ANGTOB = ANGTOB - 180.
	IF(ANGTOB.LT.-90.) ANGTOB = ANGTOB + 180.
	IF(ANGTOB.GT.90.) ANGTOB = ANGTOB - 180.
	IF(ANGTOB.LT.-90.) ANGTOB = ANGTOB + 180.
	IF(ANGTOB.LT.-90.) ANGTOB = ANGTOB + 180.
c        print*,'angtob',angtob
c        if(no_evt.eq.12115717) stop
C
C	ITEM = 'CAL_PA'
C	ok = W_ITEM_I4(ch, item, ICAL, 1, ret_size)
C	IF(ICAL.NE.0) GO TO 100
C	HAVE TO EXCLUDE CALS BY TIMING.  
c		to exclude cals
	  if(sceti4(2).ge.010000.and.sceti4(2).lt.011000) THEN
		PRINT*,'CALIBRATE'
		GO TO 100
	  endif
C
	CALL DIFFFK(XRE,YRE,ZRE,BX,BY,BZ,DDIFF,XT,YT,ZT,XI,YI,ZI,IN)
C
	TDIST = (XRE-XT)**2 + (YRE-YT)**2 + (ZRE-ZT)**2
	TDIST = SQRT(TDIST)
C
C
C	  FIND MAXIMUM IN SPECTRUM
C
	  SMAX = -1000.
	  DO N = 35,1022			! 35 IS ABOUT 2 kHz
	    IF(DBSPEC1(N).GT.SMAX) THEN
	      SMAX = DBSPEC1(N)
	      NSPMAX = N
	    ENDIF
	  ENDDO
	PRINT*,'X SPECTRUM MAX:N,F,SMAX',NSPMAX,(NSPMAX-1)*SPS/2048.,SMAX
	  FREQMAX1 = (NSPMAX-1)*SPS/2048.
C
	  SMAX = -1000.
	  DO N = 35,1022			! 35 IS ABOUT 2 kHz
	    IF(DBSPEC2(N).GT.SMAX) THEN
	      SMAX = DBSPEC2(N)
	      NSPMAX = N
	    ENDIF
	  ENDDO
	PRINT*,'Y SPECTRUM MAX:N,F,SMAX',NSPMAX,(NSPMAX-1)*SPS/2048.,SMAX
	  FREQMAX2 = (NSPMAX-1)*SPS/2048.
C
C	  AVERAGE OVER ADJACENT POINTS, THIS IS TEMPORARY FOR COMPARISON
	  FTOT = 0.
	  STOT = 0.
	  DO NN = -1,1
	    FNN = (NSPMAX-1+NN)*SPS/2048.
	    SNN = 10.**(.1*DBSPEC1(NSPMAX+NN))
	    FTOT = FTOT+FNN*SNN
	    STOT = STOT+SNN
	  ENDDO
	  FREQMAX = FTOT/STOT
	  PRINT*,'INTERPOLATED FREQ',FREQMAX
C	  WRITE(76,*) 'INTERPOLATED FREQ',FREQMAX
C
C	CRITERIA
C
	print*,'start criteria, diff=',ddiff
c	IF(DDIFF.GT.0.) GO TO 100
C	if(no_evt.ne.74375) go to 100
	call FANDBW(CH,FREQ,DBSPEC1,DBSPEC2,NPEAK,FAVR,FBW,SKEW)
	
C
	XMAX = 0.
	YMAX = 0.
	XMIN = 0.
	YMIN = 0.
	DO N = 1,2048
	  XMAX = AMAX1(XMAX,V1DATA(N))
	  YMAX = AMAX1(YMAX,V2DATA(N))
	  XMIN = AMIN1(XMIN,V1DATA(N))
	  YMIN = AMIN1(YMIN,V2DATA(N))
	ENDDO
C	print*,'xmin,max,y',xmin,xmax,ymin,ymax
	XMAX = AMAX1(XMAX,-XMIN)
	YMAX = AMAX1(YMAX,-YMIN)
C	CHANGE TO mV/m
C	 BASED ON LEFFX = 41.1, LEFFY = 3.79
	XMAX = 1000.*XMAX/41.1
	YMAX = 1000.*YMAX/3.79	
	  IFOUND = 0
C	  GO TO 100
C
	FP = 0.
	IF (DENS.GT.0.) FP = AMAX1(8.E3,9.E3*SQRT(DENS))     ! PLASMA FREQ IN KHZ
	IF(FREQMAX.LT..75*FP.OR.FREQMAX.GT.1.25*FP) THEN
	  IF(FP.NE.0.) THEN
     	   PRINT*,'FREQ NOT CLOSE TO PLASMA FREQ, FP,RATIO-',FP,FREQMAX/FP
	  ELSE
     	   PRINT*,'FREQ NOT CLOSE TO PLASMA FREQ, FP,F ',FP,FREQMAX
	  ENDIF
	  GO TO 100
	ENDIF
C
C	ADD TO BACKGROUND
C
	IBCKG = 0
	IF(IBCKG.NE.0) THEN
	  XMXLIM = .1
	  IF(XMAX.LT.XMXLIM.AND.TDS_CHANNEL.EQ.1) THEN
	    DO N = 1,1025
	      TDSFBCK(N,TDS_CHANNEL,ISPS+1) =
     1      TDSFBCK(N,TDS_CHANNEL,ISPS+1) + DBSPEC1(N)
	    ENDDO
	    NBCK(TDS_CHANNEL,ISPS+1) = NBCK(TDS_CHANNEL,ISPS+1) +1
	  ENDIF
C
	  XMXLIM = .2
	  IF(YMAX.LT.XMXLIM.AND.TDS_CHANNEL.EQ.2) THEN
	    DO N = 1,1025
	      TDSFBCK(N,TDS_CHANNEL,ISPS+1) =
     1      TDSFBCK(N,TDS_CHANNEL,ISPS+1) + DBSPEC1(N)
	    ENDDO
	    NBCK(TDS_CHANNEL,ISPS+1) = NBCK(TDS_CHANNEL,ISPS+1) +1
	  ENDIF
	ENDIF		! END OF ADD TO BACKGROUND
C
C	CALCULATE BANDWIDTH, ETC.
	NFLLIM = NSPMAX/2
	NFLLIM = MAX0(NFLLIM,NFLLIMP+1)
	NFUPLIM = (3*NSPMAX)/2
	NFUPLIM = MIN0(NFUPLIM,NFUPLIMP)
	PRINT*,'DO FREQS',FREQ(NFLLIM),FREQ(NFUPLIM)
	CONST1 = (1000./41.1)**2
	CONST2 = (1000./3.79)**2
	DO NC = 1,2
	  FCOUNT(NC) = 1.E-20
	  FAVR(NC) = 0.
	  FSTD(NC) = 0.
	  F3MOM(NC) = 0.
	  DO N = NFLLIM,NFUPLIM
	    IF(NC.EQ.1) THEN
		VLT = CONST1*10.**(.1*DBSPEC1(N))
	    ELSE
		VLT = CONST2*10.**(.1*DBSPEC2(N))
	    ENDIF
	    FCOUNT(NC) = FCOUNT(NC) + VLT
	    FAVR(NC) = FAVR(NC) + VLT*FREQ(N)
	    FSTD(NC) = FSTD(NC) + VLT*FREQ(N)**2
	    F3MOM(NC) = F3MOM(NC) + VLT*FREQ(N)**3
	  ENDDO
	  FAVR(NC) = FAVR(NC)/FCOUNT(NC)
	  FSTD(NC) = FSTD(NC)/FCOUNT(NC) - FAVR(NC)**2
	  FBW(NC) = SQRT(AMAX1(FSTD(NC),0.))
	  SKEW(NC) = F3MOM(NC)/FCOUNT(NC) - 3.*FAVR(NC)*FSTD(NC)
     1		 - FAVR(NC)**3
	print*,' bw and .02*f',fbw(nc),.02*favr(nc)
	  IF(FBW(NC).GT..02*FAVR(NC)) GO TO 100
	ENDDO
C
C	DO SAME FOR HARMONIC
C
C	CALCULATE BANDWIDTH, ETC.
	NBW = 3
	NFLLIM = 2*NSPMAX - 1 - NBW
	NFLLIM = MAX0(NFLLIM,NFLLIMP+1)
	NFLLIM = MIN0(NFLLIM,NFUPLIMP)
	NFUPLIM = 2*NSPMAX - 1 + NBW
	NFUPLIM = MIN0(NFUPLIM,NFUPLIMP)
	PRINT*,'CHECK NFUPLIM, ETC.',NFLLIM,NFUPLIM
	IF(NFUPLIM.GT.0) 
     1		PRINT*,'DO HARM FREQS',FREQ(NFLLIM),FREQ(NFUPLIM)
	DO NC = 1,2
	  F2COUNT(NC) = 1.E-20
	  F2AVR(NC) = 0.
	  F2STD(NC) = 0.
	  DO N = NFLLIM,NFUPLIM
	    IF(NC.EQ.1) THEN
		PWR = CONST1*10.**(.1*DBSPEC1(N))
	    ELSE
		PWR = CONST2*10.**(.1*DBSPEC2(N))
	    ENDIF
	    F2COUNT(NC) = F2COUNT(NC) + PWR
	    F2AVR(NC) = F2AVR(NC) + PWR*FREQ(N)
	    F2STD(NC) = F2STD(NC) + PWR*FREQ(N)**2
	  ENDDO
	  F2AVR(NC) = F2AVR(NC)/F2COUNT(NC)
	  F2STD(NC) = F2STD(NC)/F2COUNT(NC) - F2AVR(NC)**2
	  F2BW(NC) = SQRT(AMAX1(F2STD(NC),0.))
	ENDDO
C
c*************
c	CALL FANDBW(SPS,TDS_CHANNEL,DBSPEC1,DBSPEC2,NPEAK,FPK,BWPK,AMPK)	
c	print*,'fandbw',npeak,fpk(1),bwpk(1),ampk(1)
c	if(1) stop
c*************
C	FCE = .028*BMAG
C	BETA = .402*(TEMPE+TEMPI)/BMAG**2
C
c	IF(FAVR(1).GT..8*FP) THEN
	  PRINT*,'EVENT NO.',NO_EVT,' FREQ', FAVR(1),' MAX',
     1		MAXCH1,MAXCH2
c	  IFOUND = 0
c	ENDIF
C
C*****	CHECKING
c	PRINT*,'VMATRIX',VMATRIX(1,1),VMATRIX(1,2),VMATRIX(2,2)
c	PRINT*,'EVALS,AXS,PX,PY',EVAL(1),EVAL(2),AXMAJOR,AXMINOR,
c     1		PX,PY
C*****
C
C	PLASMA PARAMETERS, ETC.

C	   write(s,'(i8.8,i6.6)',iostat=ios) s_scet(1), s_scet(2)
C	   s_scet = s(1:4)//'/'//s(5:6)//'/'//s(7:8)//' '//
C	1	s(9:10)//':'//s(11:12)//':'//s(13:14)
C
c	print date, time, event no., freq, Px, Py, Vy/Vx, angles of
c		antenna,B, phiB, mag B or fce, fp
c
c	PRINT*, 'PXTOT,PYTOT,RATIO',PXTOT,PYTOT,XYRATIO
C
	print*,'write 90'
C	WRITE FREQ,AMP,BW FOR CH 1 AND 2 FUNDAMENTAL, AMP AND BW FOR
C	    HARMONIC, DIFF, DISTANCE TO TANGENT PT., FP, ANGLE TO B,
	WRITE(90,1011) sceti4,NO_EVT,EVENT(1:1),.001*FAVR(1),XMAX,
     1    .001*FBW(1),.001*FAVR(2),YMAX,.001*FBW(2),.001*FP,
     2	  ANGTOB,DDIFF,TDIST,FCOUNT(1),FCOUNT(2),F2COUNT(1),F2COUNT(2)
C
 1011	format(I10,I7,I10,A2,F6.2,E9.2,F6.3,F6.2,E9.2,F6.3,F6.2,F7.1,
     1	F6.1,F6.1,4E9.2)
C
        WRITE(80,1080) sceti4,NO_EVT,ERTDAY,EVENT(1:1),.001*FAVR(1),
     1    XRMS,.001*FBW(1),.001*FAVR(2),YRMS,.001*FBW(2),.001*FP,
     2    ANGTOB,AXRATIO,THETAB,DDIFF,XRE,YRE,ZRE
C
 1080   format(I10,I7,I10,I3,A2,F6.2,E9.2,F6.3,F6.2,E9.2,F6.3,F6.2,F7.1,
     1  F6.2,F6.1,4F7.1)
C
	IF(FAVR(1).LT.25000.) THEN
	  TPWR = FCOUNT(1) + FCOUNT(2)
	  T2PWR = F2COUNT(1) + F2COUNT(2)
	  WRITE(91,1012) sceti4,NO_EVT,EVENT(1:1),.001*FAVR(1),
     1	   .001*F2AVR(1),F2AVR(1)/FAVR(1),FCOUNT(1),F2COUNT(1),
     2	   TPWR,T2PWR
 1012	format(I10,I7,I10,A2,2F6.2,F6.3,4E11.3)
	ENDIF
C
	IF(OK.NE.82) GO TO 100
	CLOSE(UNIT=90)
	return
	end
	SUBROUTINE FANDBW(CH,FREQ,DBX,DBY,NPEAK,FAVR,FBW,SKEW)
C
	COMMON /NOISE_LIM/ XNS
	INTEGER CH
	DATA SPS /120000./
	REAL AMP(2),DBX(1025),DBY(1025),DBSUM(1025),FREQ(1025)
	REAL FCOUNT(2),FAVR(2),FSTD(2),F3MOM(2),FBW(2),SKEW(2)
	REAL X(25)
	DATA XNS /-120./
	DATA NFLLIMP,NFUPLIMP /4, 1024/
	DATA SPS /120000./
C
C	ROUTINE TO FIT FUNCTION TO LOG OF SPECTRAL POWER
C
C	CONST1 = (1./41.1)**2		! DO IN (VOLTS/M)**2
C	CONST2 = (1./3.79)**2
	LEFFX = 41.1
	LEFFY = 3.79
C
C	CALCULATE MAGNITUDE OF FIELD
C
	DO N = 1,1025
	  TPWR = 10.**(.1*DBX(N))/LEFFX**2 + 
     1		10.**(.1*DBY(N))/LEFFY**2
	  DBSUM(N) = 10.*ALOG10(TPWR)
	ENDDO
C
C	FIND FREQUENCY OF HIGHEST PEAK BETWEEN N1 AND N2
C
	N1 = 4			! LOWEST FREQ 177 HZ
	N2 = 1024
	INC = 1
C
	SMAX = DBSUM(N1)
	SSAVE = DBSUM(N1)
	NPEAK = N1
	SLIM = AMAX1(XNS,SSAVE-80.) 
	DO N = N1,N2
	    IF(DBSUM(N).GT.SMAX) THEN
	        SMAX = DBSUM(N)
	        NPEAK = N
	    ENDIF
	ENDDO
cdiag	PRINT*,'in FANDBW, MAX: N,F,SMAX',NPEAK,FREQ(NPEAK),SMAX
	FPK = FREQ(NPEAK)
	print*,'freq check',(npeak-1)*sps/2048.,fpk
C
C	LIMIT TO FAIRLY LARGE PEAKS
C
c	IF(DBSUM(NPEAK).LT.(XNS+20.)) GO TO 100
c	IF(NPK.LE.10) GO TO 100
C
C	DETERMINE DATA TO USE
C
	N2 = NPEAK+1
	SLIM = AMAX1(XNS,SMAX-25.)
	DO N = NPEAK+1,1024
	  IF(DBSUM(N).LT.SLIM) GO TO 22
	  N2 = N
	ENDDO
C
 22	N1 = NPEAK-1
	DO N = NPEAK-1,3,-1
	  IF(DBSUM(N).LT.SLIM) GO TO 24
	  N1 = N
	ENDDO
 24	CONTINUE
	print*,'fit range',n1,n2
C	if(1) stop
C
C	FIT
C
	RETURN
	END
	SUBROUTINE FANDBW2(N1,N2,FREQ,DBX,DBY,NPEAK,FAVR,FBW,SKEW)
C
	DATA SPS /120000./
	REAL AMP(2),DBX(1025),DBY(1025),DBSUM(1025),FREQ(1025)
	REAL FCOUNT(2),FAVR(2),FSTD(2),F3MOM(2),FBW(2),SKEW(2)
	DATA NFLLIMP,NFUPLIMP /4, 1024/
C
	CONST1 = (1./41.1)**2		! DO IN (VOLTS/M)**2
	CONST2 = (1./3.79)**2
C
C	CALCULATE BANDWIDTH, ETC.
C
	NFLLIM = NPEAK/2
	NFLLIM = MAX0(NFLLIM,NFLLIMP+1)
	NFUPLIM = (3*NPEAK)/2
	NFUPLIM = MIN0(NFUPLIM,NFUPLIMP)
	IF(N1.GT.N2) THEN
	   N2T = N1
	   N1T = N2
	ELSE
	   N1T = N1
	   N2T = N2
	ENDIF
C*********
C	THIS CHOICE DID NOT WORK WELL.  IT GAVE NARROW BANDWIDTHS
C	FOR SOME CASES WHEN THE BANDWIDTH WAS WIDE.  FOR EXAMPLE, IT
C	GAVE A BANDWIDTH OF 63.7 HZ FOR 19960421 # 13042040
C*********
C
	NPEAK = N1T
	VMAX = 0.
	DO NC = 1,2
	  FCOUNT(NC) = 1.E-20
	  FAVR(NC) = 0.
	  FSTD(NC) = 0.
	  F3MOM(NC) = 0.
	  DO N = N1T,N2T
	    IF(NC.EQ.1) THEN
		VLT = CONST1*10.**(.1*DBX(N))
	    ELSE
		VLT = CONST2*10.**(.1*DBY(N))
	    ENDIF
	    IF(VLT.GT.VMAX) THEN
	       NPEAK = N
	       VMAX = VLT
	    ENDIF
	    FCOUNT(NC) = FCOUNT(NC) + VLT
	    FAVR(NC) = FAVR(NC) + VLT*FREQ(N)
	    FSTD(NC) = FSTD(NC) + VLT*FREQ(N)**2
	    F3MOM(NC) = F3MOM(NC) + VLT*FREQ(N)**3
	  ENDDO
	  FN = FCOUNT(NC)
c	print*,'sk',fn,f3mom(nc)/fn,3.*fstd(nc)/fn**2,2.*(favr(nc)/fn)**2
c	  SKEW(NC) = F3MOM(NC)/FN - 3.*FSTD(NC)/FN**2  
c     1		 + 2.*(FAVR(NC)/FN)**3
	  FAVR(NC) = FAVR(NC)/FCOUNT(NC)
	  FSTD(NC) = FSTD(NC)/FCOUNT(NC) - FAVR(NC)**2
	  FBW(NC) = SQRT(AMAX1(FSTD(NC),0.))
	  SKEW(NC) = F3MOM(NC)/FCOUNT(NC) - 3.*FAVR(NC)*FSTD(NC)
     1		 - FAVR(NC)**3
	ENDDO
C
	RETURN
	END
