C+ C NAME: C TheFit C PURPOSE: C Subtracts a zodiacal light model from the raw zodiacal light data. C CATEGORY: C Data processing C CALLING SEQUENCE: program THEFIT C INPUTS: C Helios data file containing data for one color and one filter. C OUTPUTS: C Graphics to screen; binary Helios file with normalized data; ASCII file C THEFIT.DAT containing results of LSQ fits. C CALLS: C BadR4, BadR8, AskYN, AskWhat, AskI4, AskR4, AskR8, TimeSystem, LapsedTime C iSetFileSpec, iPutFileSpec, iGetFileSpec, bOpenFile, iFreeLun C iHOSInfo, HOSRead, HOSWrite, Say, DustAsymmetry, IndexR8 C LSQ_REDUCED_FIT, LSQ_CHECK_FIT, LSQ_SLEDGE_HAMMER, LSQ_UITBIJTER C iSetTermAndPlot, iSetTermOff, iSwitchGraphicsOn, iSwitchGraphicsOff, ClearWindow, C iSetSysNormal, GetCharInfo, XYLabel, iGetDeviLimits, USRDF, iSetSysUser, NOUT C iSetRightJustify, iSetBottomJustify, iSetTopJustify, iSetCenterJustify C iSetHorizJustify, iSetVertJustify, Hardcopy, GinBar, RemoveTag, bIsXevent C DrawBox, Str2StrSet, Str2Str, Int2Str, Flt2Str, Move, Draw, MoveDraw C INCLUDE: C include 'e9_hos.h' !(NOT USED ANYMORE) include 'hos_e9.h' include 'filparts.h' include 'str2str_inc.h' C COMMON BLOCKS: C common /FIT/ FLNG(MXC) ,CLNG(MXC) ,SLNG(MXC) ,FINT(MXC), NFIT C common /REJECT/ FLNGREJ(MXC),CLNGREJ(MXC),SLNGREJ(MXC),FINTREJ(MXC),NREJ C common /FITPARAM/ E,A,B,PHI,CHI, /FITPOW/ POW C common /BESTFIT/ A_IN,B_IN,PHI_IN,POW_IN,CHI_IN,DEV_IN,DSIG_IN,IPLOT C SIDE EFFECTS: C RESTRICTIONS: C PROCEDURE: C >> The zodiacal light curve to is obtained by fitting the function C I = A+B*(1+E*cos(LNG-LNGP+PHI))**POW to the data. E and LNGP are C ellipticity and ecliptic longitude of the perihelion of the HELIOS C orbit; LNG is the spacecraft ecliptic longitude. The parameters A,PHI C and POW are fitted in a least square sense. C >> The input data files contain raw data (as taken from tape by HERDA) C which are processed first using CONNECT. C >>>> USE ONLY OUTPUT FILES OF CONNECT AS INPUT FILES FOR THIS PROGRAM. C These files contain chronologically ordered data for one particular C color and one particular filter and may contain data for up to three C photometers. The names for these files have the form XyrLF_doy.ZLD C >> After substraction of the intensity of the zodiacal dust cloud, the C residual data are stored on the files XyrLF_doy.DAT C MODIFICATION HISTORY: C- parameter (MXC=2500) parameter (NS=16) ! Start sector (modified) parameter (NST=1-NS) ! End sector (modified) parameter (NSTOT=NS-NST+1) ! Total # sectors integer ICNT(0:2) /3*0/ integer NN(MXC) integer PP(MXC) integer CC(MXC) integer FF(MXC) integer LIND(MXC) ! Used for indexing integer Lst(3,11) integer Int2Str integer Str2Str integer Str2StrSet integer Flt2Str real TT(MXC) real RR(MXC) real LL(MXC) real ZZ(NST:NS,MXC) real REST(MXC) real ZCORE(NSTOT) real ELL(2) /.521757,.543814/ ! Orbit ellipticity real PHI_SC(2) /257.8477,294.2342/ ! Perihelion longitude real*8 FLNG real*8 CLNG real*8 SLNG real*8 FINT real*8 FLNGREJ real*8 CLNGREJ real*8 SLNGREJ real*8 FINTREJ real*8 RLNG real*8 E real*8 A real*8 B real*8 PHI real*8 POW real*8 CHI real*8 A_OLD real*8 B_OLD real*8 PHI_OLD real*8 POW_OLD real*8 CHI_OLD real*8 A_IN real*8 B_IN real*8 PHI_IN real*8 POW_IN real*8 CHI_IN real*8 Bad8 common /FIT/ FLNG(MXC) ,CLNG(MXC) ,SLNG(MXC) ,FINT(MXC),NFIT common /REJECT/ FLNGREJ(MXC),CLNGREJ(MXC),SLNGREJ(MXC),FINTREJ(MXC),NREJ common /FITPARAM/ E,A,B,PHI,CHI, /FITPOW/ POW common /BESTFIT/ A_IN,B_IN,PHI_IN,POW_IN,CHI_IN,DEV_IN,DSIG_IN,IPLOT logical bSym logical bBatch logical bGraphics /.TRUE./ logical bIsXEvent logical bOneCC logical bOneFF character cFileIn*80 character cFileOut*80 character cFileAsc*10 /'thefit.dat'/ character cStr*80 character cSay*6 /'TheFit'/ real*8 dsind real*8 dcosd integer time(2) bOneCC(kC) = kC .eq. HOS__C_1 .or. kC .eq. HOS__C_2 .or. kC .eq. HOS__C_3 bOneFF(kF) = kF .eq. HOS__F_1 .or. kF .eq. HOS__F_2 .or. kF .eq. HOS__F_3 .or. kF .eq. HOS__F_4 .or. kF .eq. HOS__F_5 Bad4 = BadR4() Bad8 = BadR8() call AskYN('Batch mode$no',bBatch) if (bBatch) call AskYN('Graphics output$no',bGraphics) if (bGraphics) then I = iSetTermAndPlot(-1,1) if (bBatch) I = iSetTermOff() end if if (iHOSInfo(0,cFileIn,iRecl,iSc,iYr,TMin,TMax,kPCFall,Lst) .eq. 0) & call Say(cSay,'W','StopPlay','Empty input file or read error: '//cFileIn) if (iand(kPCFall,HOS__NONORM) .eq. 0) & call Say(cSay,'W','StopPlay','Data are already normalized: '//cFileIn) kC = iand(kPCFall,HOS__C_ALL) if (.not. bOneCC(kC)) call Say(cSay,'W','StopPlay','Only one color allowed: '//cFileIn) !------- ! Read 15/30 degree data records : ICNT(1) = # 16/31 deg records ! We explicitly exclude pB data (filter 5) kPCF = kPCFall kPCF = kPCF-iand(kPCF,HOS__P_ALL)+HOS__P_12+HOS__SWAP_SECT ! Phot 1&2; swap sectors, single color kPCF = kPCF-iand(kPCF,HOS__F_ALL)+HOS__F_POL+HOS__F_CLR ! Exclude pB call HOSRead(kPCF,cFileIn,MXC,ICNT(1),NN,TT,PP,CC,FF,RR,LL,NSTOT,ZZ) kF = iand(kPCF,HOS__F_ALL) if (.not. bOneFF(kF)) call Say(cSay,'W','StopPlay','Only one filter allowed: '//cFileIn) if ((kF .eq. HOS__F_4 .or. ICNT(1) .eq. 0) .and. iand(kPCFall,HOS__P_3) .ne. 0) then !------- ! Read 90 deg data records : ICNT(2)-ICNT(1) = # 90 deg records kPCF = kPCFall kPCF = kPCF-iand(kPCF,HOS__P_ALL)+HOS__P_3 ! Phot 3, single color I = ICNT(1)+1 call HOSRead(kPCF,cFileIn,MXC,ICNT(2),NN(I),TT(I),PP(I),CC(I),FF(I),RR(I),LL(I),NSTOT,ZZ(NST,I)) end if ICNT(2) = ICNT(1)+ICNT(2) !------- ! # sectors: NSTOT for 15-30 deg and 2 for 90 deg IBEG = 1 if (ICNT(1) .eq. 0) IBEG = 2*NSTOT+1 ! No 15/30 degree photom IEND = 2*NSTOT if (ICNT(2) .ne. ICNT(1)) IEND = IEND+2 ! If 90 degree photom present E = dble(ELL(iSc)) ! Orbit ellipticity call AskYN('Transform to dust cloud plane of symmetry$yes',bSym) if (bSym) call DustAsymmetry(2,iSc,ICNT(2),PP,LL,NSTOT,ZZ) call AskWhat('Fit photometer/sector combinations: All,Single$1$1',I) if (I .eq. 1) then I = iPutFileSpec(FIL__NAME,FIL__TYPE,cFileAsc) I = iGetFileSpec(0,FIL__TYPE,cFileOut) iRecl = 0 if (bOpenFile(OPN__TEXT+OPN__NEW+OPN__TRYINPUT+OPN__ONEPASS+OPN__STOP,iU,cFileOut,iRecl)) continue c iU = iGetLun(cFileOut) c open (iU, file=cFileOut, status='NEW', carriagecontrol='LIST') write (iU,'(A)') 'FILE '//cFileIn(:itrim(cFileIn)) else call AskWhat('Photometer: 1,2,3',IPHOT) if (IPHOT .eq. 3) then I = 1 call AskWhat('Fit: Intensity, pB',I) ISECT = NST if (I .eq. 2) ISECT = NST+1 else call AskI4('Sector: 1,..,32$1$32',ISECT) ISECT = ISECT-16-sign(1,ISECT-17)*16 ! Modified sector number end if IBEG = NSTOT*(IPHOT-1)+ISECT+NS IEND = IBEG ! Just one photometer/sector combination selected end if if (bGraphics) I = iSwitchGraphicsOn() IPHOT = -1 do IPLACE=IBEG,IEND ! Loop over selected photom/sector combinations I = IPHOT ! Photometer for previous 'sector' IPHOT = 1+(IPLACE-1)/NSTOT if (I .ne. IPHOT) PHI_IN = Bad8 ! Changing photometer if (PHI_IN .ne. Bad8 .and. (abs(PHI_IN) .gt. 10 .or. & POW_IN .lt. 1 .or. POW_IN .gt. 5) ) PHI_IN = Bad8 ISECT = IPLACE-NSTOT*(IPHOT-1)-NS ! Modified sector number write (*,'(/,A,I2,A,I3)') ' Photometer ',IPHOT,'; Sector ',ISECT I2 = (IPHOT+1)/2 ! I2=1 for phot 1 & 2; I2=2 for phot 3 I1 = ICNT(I2-1)+1 I2 = ICNT(I2) NFIT = 0 PMIN0 = 1E10 PMAX0 = -1E10 J = 0 RLNG = 180.-PHI_SC(iSc) do I=I1,I2 if (I .gt. 1 .and. LL(I) .lt. LL(I-1)) J = J+1 if (PP(I) .eq. IPHOT .and. ZZ(ISECT,I) .ne. Bad4) then NFIT = NFIT+1 FLNG(NFIT) = LL(I)+RLNG ! SC longitude relative to perihelion CLNG(NFIT) = dcosd(FLNG(NFIT)) SLNG(NFIT) = dsind(FLNG(NFIT)) FLNG(NFIT) = FLNG(NFIT)+J*360. FINT(NFIT) = ZZ(ISECT,I) PMIN0 = min(ZZ(ISECT,I),PMIN0) PMAX0 = max(ZZ(ISECT,I),PMAX0) end if end do ISIG = 0 ! Automatic pilot A_IN = Bad8 ! No fit available yet NREJ = 0 if (NFIT .le. 0) go to 999 ! No points to fit call IndexR8(1,NFIT,1,MXC,FLNG,LIND) FLNGMIN = FLNG(LIND(1)) ! Minimum longitude FLNGMAX = FLNG(LIND(NFIT)) ! Maximum longitude XMIN = FLNGMIN-1 XMAX = FLNGMAX+5. ! Needed for plotting NLNG = FLNGMAX-XMIN call TimeSystem(time) call Say(cSay,'I','Init','Initialization 1:') IPLOT = -1 ! Store next file with LSQ_CHECK_FIT if (PHI_IN .ne. Bad8) then ! Try final fit for previous sector first PHI = PHI_IN POW = POW_IN call LSQ_REDUCED_FIT(POW) if (POW .ne. Bad8) call LSQ_CHECK_FIT else call Say(cSay,'I','Init','No previous fit available for initialization') end if PHI = 0.d0 POW = 2.3d0 call LSQ_REDUCED_FIT(POW) if (POW .eq. Bad8) then ! Whole sector flagged A_IN = Bad8 go to 999 end if call LSQ_CHECK_FIT A_OLD = A_IN B_OLD = B_IN PHI_OLD = PHI_IN POW_OLD = POW_IN CHI_OLD = CHI_IN DEV_OLD = DEV_IN DSIG_OLD = DSIG_IN if (IBEG .eq. IEND) then ! Only 1 phot/sect fitted ISIG = 3 ! Enter initial fit manually go to 333 end if !--------- ! We have a lift off. LSQ_CHECK_FIT stores the values for the fit ! parameters in the set of _IN variables. Intermediate results are ! regularly compared with the values in the _IN variables. If the ! intermediate fit is an improvement, its parameter values replace ! the old _IN values (see LSQ_CHECK_FIT). Thus the _IN set at all ! times represents the best available fit. A = A_IN B = B_IN PHI = PHI_IN POW = POW_IN CHI = CHI_IN TLAP = LapsedTime(time)/1000 write (*,*) 'Sledgehammer 1: ',TLAP,' (',TLAP,')' call LSQ_SLEDGEHAMMER call Say(cSay,'I','End','End of sledgehammer 1') if (POW .eq. Bad8) then call Say(cSay,'I','Init','Back to square one') A_IN = A_OLD B_IN = B_OLD PHI_IN = PHI_OLD POW_IN = POW_OLD CHI_IN = CHI_OLD DEV_IN = DEV_OLD DSIG_IN = DSIG_OLD else C TLAP1 = TLAP C TLAP = LapsedTime(time)/1000 C write (*,*) 'Fine tuning 1 : ',TLAP,' (',TLAP-TLAP1,')' C call LSQ_FIT C write (*,*) 'End of fine tuning' end if 333 if (ISIG .le. 1) then ! ISIG=0 (autopilot) and ISIG=1: Flag bad points DSIG = 3. ! Points worse than 3 sigma will be rejected if (ISIG .eq. 1) then ! Manipulate fit produced by automatic pilot DSIG = DSIG_IN call AskR4('Rejection level (sigma)',DSIG) end if I = NREJ ! Throw out the bad points worse than DSIG sigma call LSQ_UITBIJTER(E,DSIG) if (I .ne. NREJ) then ! If points have been thrown out write (*,*) NREJ,' point(s) rejected . ',NFIT,' points left.' NREJ_OLD = I A_OLD = A_IN B_OLD = B_IN PHI_OLD = PHI_IN POW_OLD = POW_IN CHI_OLD = CHI_IN DEV_OLD = DEV_IN DSIG_OLD = DSIG_IN A = A_IN B = B_IN PHI = PHI_IN POW = POW_IN IPLOT = -1 call Say(cSay,'I','Init','Initialization 2:') call LSQ_CHECK_FIT ! Update CHI TLAP1 = TLAP TLAP = LapsedTime(time)/1000 write (*,*) 'Sledgehammer 2: ',TLAP,' (',TLAP-TLAP1,')' call LSQ_SLEDGEHAMMER call Say(cSay,'I','End','End of sledgehammer 2') if (POW .eq. Bad8) then write (*,*) 'Back to square one' A_IN = A_OLD B_IN = B_OLD PHI_IN = PHI_OLD POW_IN = POW_OLD CHI_IN = CHI_OLD DEV_IN = DEV_OLD DSIG_IN = DSIG_OLD C else C TLAP1 = TLAP C TLAP = LapsedTime(time)/1000 C write (*,*) 'Fine tuning 2 : ',TLAP,' (',TLAP-TLAP1,')' C call LSQ_FIT C write (*,*) 'End of fine tuning 2' end if end if TLAP1 = TLAP TLAP = LapsedTime(time)/1000 write (*,*) 'End of fit: ',TLAP,' (',TLAP-TLAP1,')' else if (ISIG .eq. 2) then ! Backstep to previous fit do I=1,NREJ-NREJ_OLD ! Pick up last set of rejected points I1 = NFIT+I I2 = NREJ_OLD+I FLNG(I1) = FLNGREJ(I2) CLNG(I1) = CLNGREJ(I2) SLNG(I1) = SLNGREJ(I2) FINT(I1) = FINTREJ(I2) end do NFIT = NFIT+NREJ-NREJ_OLD NREJ = NREJ_OLD call IndexR8(1,NFIT,1,MXC,FLNG,LIND) FLNGMIN = FLNG(LIND(1)) FLNGMAX = FLNG(LIND(NFIT)) A_IN = A_OLD B_IN = B_OLD PHI_IN = PHI_OLD POW_IN = POW_OLD CHI_IN = CHI_OLD DEV_IN = DEV_OLD DSIG_IN = DSIG_OLD else if (ISIG .eq. 3) then ! Enter your own fit call AskR8('Phi',PHI) call AskR8('Power',POW) call LSQ_REDUCED_FIT(POW) IPLOT = -1 call LSQ_CHECK_FIT end if if (bGraphics) then PMIN = PMIN0 ! Plot results if there is something to plot PMAX = PMAX0 do I=1,NLNG REST(I) = A_IN+B_IN*(1.+E*dcosd(XMIN+I+PHI_IN))**POW_IN PMIN = min(REST(I),PMIN) PMAX = max(REST(I),PMAX) end do call ClearWindow I = iSetSysNormal() I = GetCharInfo(xCh,yCh) iStr = Str2StrSet(STR__NOTRIM) I = 0 I = I+Str2Str(cFileIn(:itrim(cFileIn)) ,cStr(I+1:)) I = I+Str2Str(' (phot=' ,cStr(I+1:)) I = I+Int2Str(IPHOT ,cStr(I+1:)) if (IPHOT .ne. 3) then I = I+Str2Str('-sect=' ,cStr(I+1:)) I = I+Int2Str(mod(ISECT+31,32)+1 ,cStr(I+1:)) end if I = I+Str2Str(')' ,cStr(I+1:)) call XYLabel(1.,100.-yCh,cStr(:I)) I = 0 I = I+Str2Str('A=',cStr(I+1:)) I = I+Flt2Str(sngl(A_IN),-2,cStr(I+1:)) I = I+Str2Str(', B=',cStr(I+1:)) I = I+Flt2Str(sngl(B_IN),-2,cStr(I+1:)) I = I+Str2Str(', PHI=',cStr(I+1:)) I = I+Flt2Str(sngl(PHI_IN),-2,cStr(I+1:)) I = I+Str2Str(', POW=',cStr(I+1:)) I = I+Flt2Str(sngl(POW_IN),-3,cStr(I+1:)) I = I+Str2Str(', DEV=',cStr(I+1:)) I = I+Flt2Str(DEV_IN,-5,cStr(I+1:)) iStr = Str2StrSet(iStr) call XYLabel(1.,100.-2.25*yCh,cStr(:I)) I = iGetDeviLimits(XL,YL) call USRDF(0.1*XL,0.15*YL,0.9*XL,0.9*YL,XMIN,PMIN,XMAX,PMAX) I = iSetSysUser() I = GetCharInfo(xCh,yCh) call MoveDraw(XMIN,PMIN,XMAX,PMIN) call MoveDraw(XMIN,PMIN,XMIN,PMAX) iH = iSetRightJustify() iV = iSetBottomJustify() call Move(XMIN,PMIN) call NOUT(nint(PMIN)) I = iSetTopJustify() call Move(XMIN,PMAX) call NOUT(nint(PMAX)) I = iSetCenterJustify() I = int(FLNGMIN)/30 if (FLNGMIN .gt. 0.) I = I+1 XD = I*30 do while (XD .lt. FLNGMAX) call MoveDraw(XD,PMIN,XD,PMIN-.8*yCh) call NOUT(nint(XD)) XD = XD+30 end do XD = .1*xCh YD = .1*yCh do I=1,NFIT ! Plot data points c call DOT(sngl(FLNG(I)),sngl(FINT(I))) XX = FLNG(I) YY = FINT(I) call MoveDraw(XX-.3*XD,YY,XX+.3*XD,YY) call MoveDraw(XX,YY-.3*YD,XX,YY+.3*YD) end do do I=1,NREJ ! Accentuate points excluded before fitting DUM1 = sngl(FLNGREJ(I)) DUM2 = sngl(FINTREJ(I)) call DrawBox(DUM1-XD,DUM2-YD,DUM1+XD,DUM2+YD) end do call Move(XMIN+1,REST(1)) ! Draw least square fit do I=2,NLNG call Draw(XMIN+I,REST(I)) end do iH = iSetHorizJustify(iH) iV = iSetVertJustify(iV) if (bBatch) then call Hardcopy else cStr = '1-Flag,2-Undo,3-Cheat,4-Hardcopy,5-OK,6-NOK' ISIG = 4 call GinBar(1,ISIG,0.,0.,100.,cStr) do while (ISIG .eq. 4 .or. bIsXEvent(ISIG)) if (ISIG .eq. 4) then call Hardcopy call RemoveTag(cStr,',','4-Hardcopy') end if call GinBar(1,ISIG,0.,0.,100.,cStr) end do end if if (ISIG .eq. 6) A_IN = Bad8 ! Reject fit: flag phot/sector if (ISIG .le. 3) go to 333 ! Repeat fit end if if (IBEG .eq. IEND) then iU = iFreeLun(-iU) ! Close and delete call Say(cSay,'S','StopPlay','Exit') c write (*,*) 'Final fit:' c if (A_IN .eq. Bad8) A_IN = 0 c write (*,'(5(A,E11.4))') ' A=',A_IN,' B=',B_IN,' PHI=',PHI_IN, c & ' POW=',POW_IN,' DEV=',DEV_IN c stop ' ' end if 999 continue if (A_IN .ne. Bad8) then if (POW_IN .lt. 1 .or. POW_IN .gt. 5) then do I=I1,I2 if (PP(I) .eq. IPHOT) ZZ(ISECT,I) = Bad4 end do write (iU,'(I1,A,I2,1X,2(A,E11.4,1X),2(A,F8.3,1X),A,F6.1,A)') & IPHOT,'-',mod(ISECT+31,32)+1, & 'A=',A_IN,'B=',B_IN,'PHI=',PHI_IN,'POW=',POW_IN,'DEV=',DEV_IN,' FLAGGED' else RLNG = 180.-PHI_SC(iSc) do I=I1,I2 ! Substract least square fit if (PP(I) .eq. IPHOT .and. ZZ(ISECT,I) .ne. Bad4) then XX = LL(I)+RLNG XX = A_IN+B_IN*(1+E*dcosd(dble(XX)+PHI_IN))**POW_IN ZZ(ISECT,I) = ZZ(ISECT,I)-XX end if end do write (iU,'(I1,A,I2,1X,2(A,E11.4,1X),2(A,F8.3,1X),A,F6.1)') & IPHOT,'-',mod(ISECT+31,32)+1, & 'A=',A_IN,'B=',B_IN,'PHI=',PHI_IN,'POW=',POW_IN,'DEV=',DEV_IN end if else if (NFIT .gt. 1) then do I=I1,I2 if (PP(I) .eq. IPHOT) ZZ(ISECT,I) = Bad4 end do write (iU,'(A,I1,A,I2,3X,A)') 'Photom ',IPHOT,' Sect ', & mod(ISECT+31,32)+1,'FIT NOT SUCCESFUL; WHOLE SECTOR FLAGGED' else write (iU,'(A,I1,A,I2,3X,A)') 'Photom ',IPHOT,' Sect ', & mod(ISECT+31,32)+1,'NO DATA POINTS' end if end if end do iU = iFreeLun(iU) if (bGraphics) I = iSwitchGraphicsOff() if (bSym) call DustAsymmetry(3,iSc,ICNT(2),PP,LL,NSTOT,ZZ) I = iSetFileSpec(cFileIn) I = iPutFileSpec(FIL__TYPE,FIL__TYPE,'.dat') I = iGetFileSpec(0,FIL__TYPE,cFileOut) kPCF = HOS__SWAP_SECT ! Swap sectors call HOSWrite(kPCF,cFileIn,MXC,ICNT(2),NN,TT,PP,CC,FF,RR,LL,NSTOT,ZZ) c iRecl = 5+NSTOT ! Record length for output file c if (ICNT(1) .eq. 0) iRecl = 6 ! Only 90 degree data c c iU = iGetLun(cFileOut) c open (iU,file=cFileOut,status='NEW',access='DIRECT',recl=iRecl) c do I=1,ICNT(1) c CORE_DATA.TIME = TT(I) c CORE_DATA.PHOTOM = PP(I) c CORE_DATA.COLOR = CC(I) c CORE_DATA.FILTER = FF(I) c CORE_DATA.DISTANCE = RR(I) c CORE_DATA.POSITION = LL(I) c do J=1,NS c ZCORE(J) = ZZ(J,I) c ZCORE(NSTOT+1-J) = ZZ(1-J,I) c end do c write (iU'I) CORE_DATA,ZCORE c end do c c do I=ICNT(1)+1,ICNT(2) c DEG90_DATA.TIME = TT(I) c DEG90_DATA.PHOTOM = PP(I) c DEG90_DATA.COLOR = CC(I) c DEG90_DATA.DISTANCE = RR(I) c DEG90_DATA.POSITION = LL(I) c DEG90_DATA.INTENS = ZZ(NST,I) c DEG90_DATA.POLARIZ = ZZ(NST+1,I) c write (iU'I) DEG90_DATA c end do c iU = iFreeLun(iU) call Say(cSay,'S','StopPlay','Exit#Output files: '//cFileAsc//' and '//cFileOut) end