;+
; NAME:
;	IPS_LSQ
; PURPOSE:
;	Calculates linear least squares parameters, giving the equation of
;	the line for the given parameters X and Y and the correlation
;	coefficient for x and y parameters.
; CATEGORY:
;	Mathematical calculations
; CALLING SEQUENCE:
;	IPS_LSQ, BMJD,EMJD,NMJD,DAYPARM,X,Y,T,XSTR,YSTR,ELSTR
; INPUTS:
;	BMJD	modified julian day on which plotting program starts
;	EMJD	modified julian day on which plotting program ends
;	NMJD	number of data files used in plotting
;	X		x-axis array
;	Y		y-axis array
;	DAYPARM	specified day of the year for which data is used in
;			correlation program
; OUTPUTS:
;	output to file:
;	BMJD, EMJD, NUMBER, DAY, X, Y
;	SIGMA_B		error in the slope for line Y = BX + A
;	SIGMA_A		error in the y-intercept
;	SIGMA_Y	error in y-axis correlation parameter
; CALLS:
;	CORRELATE, TOTAL, SQRT
; MODIFICATION HISTORY:
;	August 1992, Susan Rappoport and Kari Winfield (UCSD)
;-

pro IPS_LSQ, BMJD,EMJD,NMJD,DAYPARM,X,Y,T,XSTR,YSTR,ELSTR

file = 'IPS'+strcompress(dayparm,/rem)+'.LSQ'
message, /info, 'Writing LSQ info to file '+file
openw, /get_lun, out, file

printf, out, 'MJD RANGE: ',BMJD,'-',EMJD
printf, out, 'NUMBER OF FILES READ: ',NMJD
printf, out, 'DAY: ',DAYPARM
printf, out, 'ELONG: ',ELSTR
printf, out, 'X: ',XSTR,' Y: ',YSTR,'  READ FROM MJD:'
printf, out, ' '
printf, out, 'X:','Y:','TIME (MJD)','READ FROM MJD',format='(4A15)'
printf, out, rotate( [[X],[Y],[T],[T-DAYPARM]], 4), format='(2F15.4,2I15)'

N = n_elements(X)
if N le 2 then	$
    printf, out, 'Insufficient number of data points for valid calculations' $

else begin
    Sxx = total(X*X)
    Sy  = total(Y)
    Sx  = total(X)
    Sxy = total(X*Y)

    D = N*Sxx-Sx*Sx
    A = (Sxx*Sy-Sx*Sxy)/D	; y-intercept
    B = (N *Sxy-Sx*Sy )/D	; slope

    SIGMA   = Y-(A+B*X)
    SIGMA   = total( SIGMA*SIGMA )/(N-2)
    SIGMA_Y = sqrt ( SIGMA )
    SIGMA_A = sqrt ( SIGMA*Sxx/D )
    SIGMA_B = sqrt ( N*SIGMA/D )

    printf, out, 'NUMBER OF POINTS IN ARRAY: ',N
    printf, out, 'SLOPE      : ',B,'     ERROR: ',SIGMA_B
    printf, out, 'CENTER     : (',Sx/N,',',Sy/n,')'
    printf, out, 'Y-INTERCEPT: ',A,'     ERROR: ',SIGMA_A
    printf, out, 'LEAST SQUARES LINE:  Y = ',B,'*X + ',A
    printf, out, 'UNCERTAINTY IN Y VALUES: ',SIGMA_Y
    printf, out, 'CORRELATION COEFFICIENT: ',CORRELATE(X,Y)

    A = normalfit(X,Y)
    printf, out, ' Normal fit:'
    printf, out, 'SLOPE      : ',-A(0,0)/A(2,0)
    printf, out, 'CENTER     : (',A(1,0),',',A(3,0),')'

endelse

free_lun, out

return  &  end