PRO PMTRAS_DBASE_CORCO, dbsolution, coeff
;
; For each 64-second data point in a database solution, calculates the correlation coefficient between
;   the observed and predicted intensities of identified stars.
; This should provide independent confirmation of correct star identification.
;
; dbsolution    = npt-vector of master database structures, where npt is the number of data points
; coeff         = corresponding vector of correlation coefficients
;
; 2007-Jun-08   Initial version adapted from Shane Frewen's gencorco routine.
;

;create an object made of the roll database
roll_db = hsi_roll_db_full()

;create a structure using the object roll_db with a date, here during may 15 of 2003
date = '2002-Aug-5 20:32'
interval = [0, 14*24*60*60]
roll_db_date = roll_db -> getdata(obs_time_interval = anytim(date) + interval)

;create 2D arrays identifying the observed intensity and Harvard reference Number
;for stars at a the 64 second intervals spanning our observing period.
;Also, creating an array of the spacecraft time at each interval, for comparision
;with the correlation coefficient
obs_intensity = dbsolution.stars.intensity          ; 10 x npt array
id = dbsolution.stars.id
sctime = dbsolution.sctime
;
; read the star catalog into a starcat structure.
starcat = hsi_pmtras_rd_starcat()

;make arrays composed of the useful values in the star catalogue, here the predicted
;intensity and the corresponding Harvard Reference Number
pred_intensity = starcat.predcounts
catid = starcat.hrn

;the correlation coeffecients will all be put in one array, and that array must be
;big enough to hold them, so the array is sized equal to the number of 64 second
;intervals, here labeled m
npt = N_ELEMENTS(dbsolution)
coeff = FLTARR(npt) - 2         ; preset to -2 as a flag for missing datapoints.
nstar = dbsolution.starcount
starindex = INTARR(10)
;
; Begin loop over time intervals.
;run through each interval, outputting the important information, correlation coefficient
;and corresponding arrays for both predicted and observed intensity
for j = 0, npt-1 do begin
   ;creating two arrays with the information from the the jth 64 second interval,
   ;and eliminating the zeroes within them, repeating for all intervals. Called
   ;'row' because they are one row from the 2D matrix that they were a part of
   ;If no stars are recorded, the interval is ignored.
    IF nstar[j] LT 2 OR dbsolution[j].roll_quality LT 200 THEN CONTINUE           ; Skip calculation for this point if fewer than 2 stars
    FOR n=0, nstar[j]-1 DO BEGIN
        i = WHERE(catid EQ id[n,j])
        IF N_ELEMENTS(i) NE 1 OR i LT 0 THEN MESSAGE, 'Unknown star(s)'
        starindex[n] = i
    ENDFOR
    coeff[j]  = CORRELATE(obs_intensity[0:nstar[j]-1, j], pred_intensity[starindex[0:nstar[j]-1]])
ENDFOR
RETURN
END