SUBROUTINE SAX_ACC_B1S1_Y(INCREMENT_ROUTINE) C C---------------------------------------------------------------------- C C.IDENTIFICATION: Subroutine SAX_ACC_B1S1_Y C.LIBRARY: FOTLIB C.AUTHOR: L.Chiappetti - IFCTR Milano C.VERSIONS: 0.0 - 10 May 93 - Original version C 0.1 - 23 Mar 94 - support to concatenation and break C 0.2 - 06 Feb 96 - disable EOD check, added packet counter C 0.3 - 07 May 96 - support to time windows and single recycle C 0.4 - 14 May 96 - support to instrument corrections C 0.5 - 19 Jun 96 - protect against "no data" observations C 0.6 - 08 Jul 96 - instrument from GET_GLOBAL_DEFAULT (MG/SDC) C 0.7 - 13 Aug 96 - added calculation of exposure time (ddf) C 0.8 - 12 Nov 96 - times in double precision C provisional x.x - 04 Dec 96 - PROVISIONAL patch for FOT time anomaly !!!!!! C indicated with comment !ANOMALY C 0.9 - 05 Dec 96 - handle "original" time resolutions C 0.10- 18 Mar 97 - protect against duplicated packets in C FOTs (ddf) C 0.11- 11 Jun 97 - two bugs in SAX_ACC_B1S1_Y and in C EXPOSURE_B1S1 corrected (ddf) C 0.12- 04 Aug 97 - FINAL patch for FOT time anomaly (i.e. s/c C anomaly 80) recovery C 0.13- 05 Aug 97 - inclusion of exposure calculation for C time bins (both for time series and C for time analysis programs) (ddf) C 0.14- 12 Jan 98 - removed IFRACTION debug statements C.PURPOSE: Laben-style direct mode data wrapper for main accumulation loop C.METHOD: loop on packets and events and call INCREMENT_ROUTINE C.SYNTAX: CALL SAX_ACC_B1S1_Y(INCREMENT_ROUTINE) C.PARAMETERS: EXTERNAL increment_routine : real accumulation routine C.RESTRICTIONS: this processes FOT data (fields are in bytes) C.NOTES: Direct mode data correposnd to "basic type" 1 (bt=1) C Secondary type 1 is the Laben format (packet is made by C one header and a number of events) C INCREMENT_ROUTINE is an external routine usually found C with the caller program, and is the work-horse of the C accumulation (is specific of the kind of accumulation) C.FILES: C.REFERENCES C C---------------------------------------------------------------------- C C+++++ C accumulation for Laben Fot Direct data (bt#1:st#1:du=bytes) C+++++ INTEGER ITYPE,ID,JE,I,J1,J2,IHL,IVS,IVE integer ip1,ip2,iev !ANOMALY logical reverse,previous !ANOMALY INTEGER PACKET_POINTER,EVENT_POINTER,IERR,NEV,NEVENT INTEGER LUTMP,IOBS * INTEGER LUTMP,IOBS,IFRACTION * REAL FRACTION INTEGER ISTIME,ITIMSIZ,ICOUNTER,LSBORIG,BITORIG INTEGER CURRENT_COUNTER,PREVIOUS_COUNTER DOUBLE PRECISION PREVIOUS_TIMEHI,UDOUBLE DOUBLE PRECISION PACKET_START_TIME,PACKET_END_TIME DOUBLE PRECISION TIMED, OLDTIME, TCYCLE,THR,JUMP LOGICAL Z_BREAK, DOUBLEEVENT, REJECT CHARACTER DUMMY LOGICAL FOUND,ENDOFDATA,CONVERT LOGICAL ENDWINDOW,ENDACCUM,NEWOBS,FIRST C this is an internal buffer to hold the largest "telemetry" packet buffer C--> provisionally size to 3000 bytes; could use dynamic allocation ? CHARACTER RECORD*3000 CHARACTER CVALID*2 INTEGER*2 NVALID EQUIVALENCE (NVALID,CVALID) C C equivalences to hold a single subfield of one event CHARACTER CX*4,C2(2)*2,C1(4)*1 INTEGER*4 I4 INTEGER*2 I2 EQUIVALENCE (CX,C2,C1),(CX,I4),(C2(1),I2) EXTERNAL INCREMENT_ROUTINE INCLUDE 'hcommon.inc' INCLUDE 'bincommon.inc' INCLUDE 'accumcommon.inc' INCLUDE 'syscommon.inc' INCLUDE 'timecommon.inc' INCLUDE 'wcommon.inc' C C Time cycling constant definition C TCYCLE=2.D0**32*TIMECOMMON_CONVERSION OLDTIME=0.D0 ACCUMCOMMON_TRECYCL=0.D0 C original time resolution before reformatting C if not in packetcap assume is the same CALL PKTCAP_LOOKUP('or',ITYPE,LSBORIG,DUMMY,FOUND) IF(.NOT.FOUND)THEN THR=-TIMECOMMON_CONVERSION lsborig=-16 ELSE THR=-(2.D0**LSBORIG)*TIMECOMMON_SCALE ENDIF C original time bit width before reformatting CALL PKTCAP_LOOKUP('ow',ITYPE,bitORIG,DUMMY,FOUND) IF(.NOT.FOUND)THEN WRITE(*,*)'Missing "ow" field in packetcap' C--> provisionally stop CALL Z_EXIT(999) ENDIF C jump wrongly introduced by Telespazio for s/c anomaly 80 C n.b. : lsborig is negative JUMP=(2.D0**(BITORIG+LSBORIG)) C bad time event can be rejected or corrected CALL GET_GLOBAL_DEFAULT('BADTIMES','Reject',DUMMY) CALL UPCASE(DUMMY) REJECT=DUMMY.NE.'C' IF(.NOT.REJECT)THEN WRITE(*,*)' eventual bad times will be subtracted ',JUMP,' s' ENDIF C C get some basic packet info : header length, no.of events in packet C--> all these should check if FOUND CALL PKTCAP_LOOKUP('hl',ITYPE,IHL,DUMMY,FOUND) CALL PKTCAP_LOOKUP('ni',ITYPE,NEVENT,DUMMY,FOUND) CALL PKTCAP_LOOKUP('ve',ITYPE,ID,DUMMY,FOUND) IF(FOUND)THEN C start and end byte of valid event counter CALL PKTCAP_LOOKUP('vl',ITYPE,IVS,DUMMY,FOUND) CALL PKTCAP_LOOKUP('vs',ITYPE,ID,DUMMY,FOUND) IVE=IVS+ID-1 C ELSE C--> if no valid event counter, shall derive info from invalid pointer .... ENDIF C C The location of the LSbyte of the packet counter C (to be inserted as a packetcap field?) C ICOUNTER=3 C C The location of the start time and its size in the packet header C (all other 3 times needed to calculate exposure time are assumed C to be after this field and with the same size in bytes - THIS IS C SAX SPECIFIC!) C CALL PKTCAP_LOOKUP('tl',ITYPE,ISTIME,DUMMY,FOUND) CALL PKTCAP_LOOKUP('ts',ITYPE,ITIMSIZ,DUMMY,FOUND) C C set time counters to zero C TIMECOMMON_DEAD=0.0d0 TIMECOMMON_LIVE=0.0d0 C C These are needed by exposure_b1s1 C FIRST=.TRUE. ENDWINDOW=.FALSE. NEWOBS=.FALSE. ENDACCUM=.FALSE. C set record counter to zero PACKET_POINTER =0 ACCUMCOMMON_TOTPKT =0 C--> telemetry file is opened in main .... ENDOFDATA=.FALSE. C--> is byte swap conversion needed (same for I2 and I4 provisionally) CONVERT=SYSCOMMON_ICONV.EQ.1 IF(CONVERT) . WRITE(*,*)' packet data will be read and byte-swapped ', . 'as needed' C C Needed in EXPOSURE_B1S1 C PREVIOUS_TIMEHI=TIMEHI * IFRACTION=ACCUMCOMMON_NREC/10 C C loop on packets C 1 CONTINUE * timelo=0.0d0 !ANOMALY!******* * timehi=1.D39 !ANOMALY!******* previous=.true. !ANOMALY IF (Z_BREAK())GOTO 99 PACKET_POINTER=PACKET_POINTER+1 * IF(MOD(PACKET_POINTER,IFRACTION).EQ.0) * * WRITE(*,*)(PACKET_POINTER/IFRACTION)*10,' % completed' C this means it is the last packet IF(PACKET_POINTER.GT.ACCUMCOMMON_NREC)GOTO17 C--> error policy to be implemented READ(ACCUMCOMMON_LU,REC=PACKET_POINTER,IOSTAT=IERR) . RECORD(1:ACCUMCOMMON_RECL) C C Here we check if the end time of the packet is inside C the accumulation time window (if not we read the next packet C with noop) - assuming that it is at ISTIME+ITIMSIZ C CX=RECORD(ISTIME+ITIMSIZ:ISTIME+2*ITIMSIZ-1) IF(CONVERT)CALL SWAPI4(I4,1) PACKET_END_TIME=UDOUBLE(I4)*TIMECOMMON_CONVERSION+ . ACCUMCOMMON_TRECYCL TIMED=PACKET_END_TIME ip2=i4 !ANOMALY C C Another FOT anomaly: duplicated packets C Patch to filter duplicate packets using the packet counter C (only 8 LSB) C C1(1)=RECORD(ICOUNTER:ICOUNTER) CURRENT_COUNTER=ICHAR(C1(1)) IF(CURRENT_COUNTER.EQ.PREVIOUS_COUNTER) . GOTO 1 PREVIOUS_COUNTER=CURRENT_COUNTER C C A time cycling can be anywhere.... C IF(PACKET_END_TIME.LT.OLDTIME)THEN TIMED=PACKET_END_TIME+TCYCLE C try saving the time as well in correct units ! PACKET_END_TIME=PACKET_END_TIME+TCYCLE ENDIF C C (A) check if packet is entirely before start of current window IF(TIMED.LT.TIMELO)THEN C C If there is a time cycling, we update TRECYCL C (ONLY if we loop back to label 1) C IF(PACKET_END_TIME.LT.OLDTIME) . ACCUMCOMMON_TRECYCL=ACCUMCOMMON_TRECYCL+TCYCLE OLDTIME=TIMED GO TO 1 ENDIF Cddf C Check if the start time of the new packet is inside the C current time window. If not, advance the window pointer C suitably. C CX=RECORD(ISTIME:ISTIME+ITIMSIZ-1) IF(CONVERT)CALL SWAPI4(I4,1) PACKET_START_TIME=UDOUBLE(I4)*TIMECOMMON_CONVERSION+ . ACCUMCOMMON_TRECYCL ip1=i4 !ANOMALY Cddf reverse=ip2.lt.ip1 !ANOMALY * write(*,1111)PACKET_POINTER,ip1,ip2,reverse !ANOMALY 1111 format(i7,1x,z8,1x,z8,1X,l) !ANOMALY C C A detailed description of the logic of this block-if C is included below in the code C IF(PACKET_START_TIME.GT.TIMEHI)THEN 18 CONTINUE WCOMMON_WPTR=WCOMMON_WPTR+1 ENDWINDOW=.TRUE. IF(WCOMMON_WPTR.LE.WCOMMON_N)THEN TIMELO=WCOMMON_START(WCOMMON_WPTR+WCOMMON_SOFF) TIMEHI=WCOMMON_END(WCOMMON_WPTR+WCOMMON_EOFF) IF(PACKET_START_TIME.GT.TIMEHI)THEN ENDWINDOW=.FALSE. NEWOBS=.TRUE. GO TO 18 ENDIF ELSE GO TO 99 ENDIF Cddf C Store the new window_timehi for use in exposure_b1s1 C (only if we are using the first good time window). C Loop back if the packet end time is before the current time C window C IF(FIRST)THEN PREVIOUS_TIMEHI=TIMEHI FIRST=.FALSE. ENDIF IF(PACKET_END_TIME.LT.TIMELO)GO TO 1 ENDIF C IF(FIRST)FIRST=.FALSE. ACCUMCOMMON_TOTPKT=ACCUMCOMMON_TOTPKT+1 EVENT_POINTER=IHL C shall also extract from packet header the Valid Event Counter C--> according to precached information (if not :ve: build it ???) C--> this sets NVALID assuming it is an I*2 .... CVALID=RECORD(IVS:IVE) IF(CONVERT)CALL SWAPI2(NVALID,1) Cddf C C Add the calculation of the exposure time C (THIS IS SAX SPECIFIC! To be handled via global?) C CALL EXPOSURE_B1S1(ICOUNTER,ISTIME,ITIMSIZ,RECORD,CONVERT, * PREVIOUS_TIMEHI,ENDWINDOW,ENDACCUM, * NEWOBS,TCYCLE) IF(ENDWINDOW)THEN C C update upper boundary for exposure time calculation C CAUTION - There could be problems with time cycling... C ENDWINDOW=.FALSE. PREVIOUS_TIMEHI=TIMEHI ENDIF C C loop on events C DO 10 NEV=1,NVALID C C take one event in a character buffer C C extract the relevant fields C--> extraction might also imply byte swap (on Sun, HP etc. according i2,i4 packetcap info) DO 15 I=1,ACCUMCOMMON_NF J1=EVENT_POINTER+ACCUMCOMMON_START(I) J2=EVENT_POINTER+ACCUMCOMMON_END(I) CX=RECORD(j1:j2) IF(ACCUMCOMMON_SIZE(I).EQ.1)THEN C--> to be verified ACCUMCOMMON_J(I)=ICHAR(C1(1)) ELSEIF(ACCUMCOMMON_SIZE(I).EQ.2)THEN IF(CONVERT)CALL SWAPI2(I2,1) ACCUMCOMMON_J(I)=I2 ELSEIF(ACCUMCOMMON_SIZE(I).EQ.4)THEN IF(CONVERT)CALL SWAPI4(I4,1) ACCUMCOMMON_J(I)=I4 C handle time in here C--> also if an event is outside TIME range it is rejected soon and C--> is not processed for instrument corrections !! OK? IF (I.EQ.ACCUMCOMMON_TINDEX)THEN C !ANOMALY code moved below to handle s/c anomaly 80 C OBT time in current units (s, ms, us) ACCUMCOMMON_T=UDOUBLE(ACCUMCOMMON_J(I))* . TIMECOMMON_CONVERSION+ACCUMCOMMON_TRECYCL C C handling OBT recycling C * IF(ACCUMCOMMON_T.LT.OLDTIME)THEN now replaced such that C the decrease shall be larger than the original time resolution IF(ACCUMCOMMON_T-OLDTIME.LT.THR)THEN ACCUMCOMMON_T=ACCUMCOMMON_T+TCYCLE ACCUMCOMMON_TRECYCL=ACCUMCOMMON_TRECYCL+TCYCLE ENDIF C spacecraft anomaly 80 is now handled here iev=i4 C anomaly 80 is the inversion of times between two events C IN a packet, must no be tested for the FIRST event IF(PACKET_END_TIME-ACCUMCOMMON_T.LT.THR)THEN * if(previous)write(*,1113)iobs,packet_pointer,ip1,ip2, * . nev-1,iev 1112 format(I2,I5,I7,1X,z8,F20.10,' time out of range',$) 1113 format(I2,I5,' start end ',Z8,1X,Z8/7X,I7,1X,z8) write(*,1112)iobs,packet_pointer,nev,i4,ACCUMCOMMON_T C in all cases reject the first anomalous event ? C it is out of time sequence ! C in all cases correct, but possibly do not process ACCUMCOMMON_T=ACCUMCOMMON_T-JUMP C the blank writes below are necessary to advance one line C if no "corrected" message is issue if(previous)then write(*,*) previous=.false. goto 9 ENDIF IF(REJECT)THEN write(*,*) goto 9 ENDIF write(*,*)' corrected ',ACCUMCOMMON_T ENDIF C this statement has to be at this point OLDTIME=ACCUMCOMMON_T C C verify time is in wished range IF (ACCUMCOMMON_T.LT.TIMELO) GOTO 9 IF (ACCUMCOMMON_T.GT.TIMEHI) GOTO 9 ENDIF ENDIF C C test PHA is in range C test if it is in range for the OTHER parameters C--> disable range checking here and move it after correction C--> otherwise one could save time rejecting here on non-correctable quantities C--> but that will complicate the IFs and probably not gain speed * JE=ACCUMCOMMON_J(I) * IF(JE.LT.ACCUMCOMMON_IEL(I).OR. * . JE.GT.ACCUMCOMMON_IEH(I))GOTO 9 15 CONTINUE C C if requested perform correction on entire event IF(ACCUMCOMMON_CORRECT)CALL CORRECT C C then test if corrected event is in range DO 16 I=1,ACCUMCOMMON_NF C bypass time to prevent side effects of time windows IF (I.EQ.ACCUMCOMMON_TINDEX)GOTO 16 JE=ACCUMCOMMON_J(I) IF(JE.LT.ACCUMCOMMON_IEL(I).OR. . JE.GT.ACCUMCOMMON_IEH(I))GOTO 9 16 CONTINUE C C now process accepted events C--> optional counters C--> NTAKEN=NTAKEN+1 CALL INCREMENT_ROUTINE 9 CONTINUE EVENT_POINTER=EVENT_POINTER+ACCUMCOMMON_EVENTSIZE C C in case of time windows verify advance to next window needed * IF(WINDOWS)THEN IF (ACCUMCOMMON_T.GT.TIMEHI) THEN C C Set a flag for exposure time calculation C ENDWINDOW=.TRUE. Cddf C C We jump here if we need to skip more than one time C window C 19 CONTINUE WCOMMON_WPTR=WCOMMON_WPTR+1 * write(*,*)' window advance ',WCOMMON_WPTR,PACKET_POINTER C C This flag must be set outside. It can be the last C window. C IF(WCOMMON_WPTR.LE.WCOMMON_N)THEN C if not after last time window C set time range to next window TIMELO= . WCOMMON_START(WCOMMON_WPTR+WCOMMON_SOFF) TIMEHI= . WCOMMON_END(WCOMMON_WPTR+WCOMMON_EOFF) Cddf C C Iterate if more than one time window is skipped (e.g. C if data are missing) and handle this case as an C "end of observation" for the exposure time calculation. C IF(ACCUMCOMMON_T.GT.TIMEHI)THEN ENDWINDOW=.FALSE. NEWOBS=.TRUE. GO TO 19 ENDIF ELSE C if after last time window terminate * WRITE(*,*)' terminating after last time window' GOTO 99 ENDIF ENDIF C C HP-GSPC specific C IF (NEV.EQ.NEVENT-1) THEN IF (ACCUMCOMMON_INSTR.EQ.'hpgs') THEN DOUBLEEVENT=.TRUE. C--> TBI check on doubled events IF (DOUBLEEVENT.EQ..TRUE.) NEVENT=NEVENT-1 ENDIF ENDIF C 10 CONTINUE C IF event is last of last packet of last observation C set end-of-data condition C it is NO LONGER necessary to consider partially full packets as end-of-data IF(NVALID.LT.NEVENT)THEN WRITE(*,*)' Partially full packet at no. ',PACKET_POINTER * WRITE(*,*)' Valid events = ',NVALID,' Events = ',NEVENT ENDIF C C We jump here from the reading loop - just in case a time window C starts in the next observation C 17 CONTINUE C this means it is the last packet IF(PACKET_POINTER.GE.ACCUMCOMMON_NREC)ENDOFDATA=.TRUE. IF(.NOT.ENDOFDATA) GOTO 1 C C if current observation is over check if a new one in chain exists C--> change of packet type not supported C CLOSE(ACCUMCOMMON_LU) LUTMP=ACCUMCOMMON_LU 11 CONTINUE CALL GET_OBS_CHAIN(IOBS,'Keep') IF(IOBS.EQ.0)GOTO 99 CALL SAX_ACC_OPEN_TLM(LUTMP,ACCUMCOMMON_PACKET) C--> if nrec is 0 file has not been opened since it does not exist IF(ACCUMCOMMON_NREC.EQ.0)GOTO 11 C reset record counter to zero PACKET_POINTER =0 ENDOFDATA=.FALSE. * IFRACTION=ACCUMCOMMON_NREC/10 C C Tell to the exposure time subroutine that a new observation was C opened C NEWOBS=.TRUE. GOTO 1 99 continue C C Before exiting, we update the exposure time C ENDACCUM=.TRUE. CALL EXPOSURE_B1S1(ICOUNTER,ISTIME,ITIMSIZ,RECORD, * CONVERT,PREVIOUS_TIMEHI,ENDWINDOW,ENDACCUM,NEWOBS,TCYCLE) RETURN END