******************************************************************************** Subroutine xtrnmsr( x, y, z, t, hdr, pl) ******************************************************************************** C C F.MSR : MSR measures areas, velocities and widths of C galaxian profiles. C It does so in up to 5 different ways, using C three different approaches: C 1 ... width at fraction f of mean flux over (PR2,PR3) C 2 ... width at fraction f of peak flux over (PR2,PR3) C 3 ... width at fraction f of each of two horns, C searched from the inside out. C 4.... width at a fraction f of each horn, located C by simple linear interpolation between the C points which bracket the point with flux = C f times the peak, the velocity uncertainity C is estimated by fitting a polynomial between C the f1 and f2 points on each side of the C profile C 5 ... width at a fraction f of each horn, searched C by fitting a polynomial between fractions C f1 and f2 of the respective horn to each C side of the profile C C C Input : $MSR PR1 PR2 PR3 C C where PR1 = sequence of up to 5 integers 1,2,3 or 4 identifying the C mode of computation of width as listed above, in up to C 5 separate calls for measuring; e.g. if one wants to C compute (V,W) first at level f1 of mean, then at level C f3 of peak, then level f4 of peak, then level f5 of C average of two horns, then at level f of each horn C located by fitting a polynomial between levels f_1 C and f_2 of the respective peak to each side of the C profile,set PR1 to 12234; values f1,f2,f3,f4,f,f_1,f_2 C are read by MSR from the PL(1,5) elements. PL(5) has to C be coded as 0.ff_1f_2 (i.e each fraction can be specif- C ied to one digit accuracy only. so if you want to fit C between 0.8 and 0.2 of the horn and the find the 0.5 C level by interpolation, set PL(5) to 0.825. Methods 4 C and 5 will return the fitted polynomials to the Z (left C half of the profile) and T (right half of the profile) C arrays. Methods 4 and 5 can also be run interactively, C to permit the user to choose the peaks, fitting levels C etc. For interactive mode, set the appropriate pl elem- C ent < 0. i.e if you want to use measuring methods 1,2, C 4,3,2 and want to run method 4 interactively, set pl(3) C < 0.0. C C C PR2, PR3 = boundaries of interval over which area and Peak search C is carried on. C C If PR1-3 are not specified, user is prompted to enter values. The C AREA returned is the one between PR2 and PR3; however, MSR also C searches for the first occurrence of nulls from the horns out, and C it computes the AREA between those two channels; that value is C returned in PL(6), and printed as AREA if approach 3 of width C measurement is invoked. It is assumed that base.f has been run C and the RMS of the profile is in STO/RCL 9 ,element 128 and that C the smoothing applied is specified in STO/RCL 9, element 121. C C C Output: Computed parameters are displayed on the screen, C On sto/rcl 9: 139 : maximum within feature C 140 : flux integral C 141-155: velocity,width, measuring code C 156 : velocity error (methods 4&5 only) c 241-250: IFRL, left ch nr for FR level (methods 4,5) c " : IFRR, right ch nr for FR level (methods 4,5) c 251 : slope of fit at IFRL (method 5) c 252 : slope of fit at IFRR (method 5) C In Pl array: 11&12 : channel numbers at which method 1 makes C the velocity measurement. C 13&14 : channel numbers at which method 2 makes C the velocity measurement. C 15&16 : channel numbers at which method 3 makes C the velocity measurement. C 17&18 : channel numbers of the peaks found by C method 4. C 19&20 : channel numbers at which method 4 makes C the velocity measurement. C 21&22 : channel numbers of the peaks found by C method 5. C 23&24 : channel numbers at which method 5 makes C the velocity measurement. C 31&32 : the channel numbers between which the C polynomial was fit by method 4.(left C half of the profile). C 33&34 : the channel numbers between which the C polynomial was fit by method 4.(right C half of the profile). C 35&36 : the channel numbers between which the C polynomial was fit by method 5.(left C half of the profile). C 37&38 : the channel numbers between which the C polynomial was fit by method 5.(right C half of the profile). C Z&T arrays 1-Isize: fitted polynomials for left and right C side of the profiles (methods 4&5, C only). C C RG/19Dec85 C JLRB/Sep89 Sun's version. C JNC/Sep90 polynomial fitting subroutines, modification of C the area measuring algorithms C JNC/Aug91 output to Pl arrays and Z arrays, modifcation of the C polynomial routines to run interactively. Slight C fiddling with the peak search algorithm. c RG/Dec92 mods to interaction mode and constrains for max error c in velocity measure entered in areaPE and areaPF; c parms written on sto/rcl 9, 241-244 C C ******************************************************************* Include '../../anzxtrna.inc' Include '../../anzxtrnb.inc' Include '../../../src/include/ttchar.inc' C Double Precision F(5),YValue,SMx,Smax,S,Dex,f Double Precision SDV(5),WRD(5),W(5),V(5) Double Precision Wd, Fint, Chan, RMS,VRMS,SFAC Double Precision slope_right, slope_left Integer*4 I, ISy, J, Imax, il, ir, i19 Integer*4 System, K(5), jl(5), jr(5),PRTLEV Integer*4 LastC, iostatus, n, m,PRTLEV Logical answer, ask CHARACTER*8 WORD(5), DpChar, capitalize CHARACTER RPLY Character*72 string EQUIVALENCE (WRD,WORD) C CHARACTER KEY*1, Ans*1 DATA K/5*0/ VRMS =0.0 IF (capitalize(DpChar(PR1)).EQ.'HELP') THEN i = System('more ANZ_BASE/external/Arecibo/Galpac/msr.1') RETURN END IF C C Zero y-array for replotting zero-line ISY = Isize Do 10 I=1,ISY Y(I) = 0.0 10 Continue ANS = "1" C C Get the master out there to settle on choice of parms C If on Graphics teminal, default to crosshair mode for bounds C IF (PR1.EQ.0.) THEN WRITE (6,'(" Enter string of modes for measure")') READ (5,*,iostat=iostatus) PR1 if (iostatus.ne.0) then call perror(' msr> ') return endif END IF C C If terminal is a Graphics terminal, use crosshair to get bounds, else C search for nulls. C 31 If (graphics) then C OK, it is a Graphics terminal... Call graphmode(.false.) call putln(' Now flag LEFT side chan for area bounds. ') call putln(' Hit any key.') Call putln(' Flag RIGHT side channel.') call putln(' Hit any key.') call graphmode(.true.) C CALL ReadCursor(KEY,CHAN,YVALUE) CALL FlagPoint(Chan) PR2 = CHAN C move to the right... chan = chan + 10.0 CALL ReadCursor(KEY,CHAN,YVALUE) CALL FlagPoint(Chan) PR3 = CHAN call graphmode(.false.) Else C So, it is not a Graphics terminal... CALL NULLS(X,PR2,PR3) CALL FlagPoint(PR2) CALL FlagPoint(PR3) End If C IF (PR2.GE.PR3) THEN call putln(" Naah... PR2 > PR3, try again.") GO TO 31 END IF C answer = ask('Happy with the choice of bounds? ', . 'Yes or No, please.', .true.) IF (.not.answer) THEN if (graphics) then call putln( ' ' ) call putln( ' Try again...' ) call putln( ' ' ) go to 31 else WRITE (6,'(" Former choice was",2F6.0)') PR2,PR3 WRITE (6,*)" Enter new channel nrs. for bounds" READ (5,*,iostat=iostatus) PR2,PR3 if (iostatus.ne.0) then call perror(' msr> ') return endif endif END IF C C Go to get the max flux C SMAX = SMx( X, PR2, PR3 ) C C Now decode PR1 and decide which routines to call C S = PR1 J = 1 DO 14 I=1,5 DEX = 10.**(5-I) K(J) = S/DEX F(J) = PL(I) S = S - Dble(K(J))*DEX IF (S.LT.0.5) GO TO 12 IF (K(J).EQ.0) GO TO 14 J = J + 1 14 CONTINUE 12 IMAX = J C Obtain RMS from STO/RCL 9, element 128 CALL Recall(Y,Isize,9) RMS = Y(128) C Compute the factor by which the smoothing has degraded the C velocity resolution, N is the box smoothing length, M is C number of Hanning smooths J = IDINT(Y(121)) N = J/10 M = J - N*10 SFAC = SQRT( N**2 + M*(1.56**2)) IF (SFAC.EQ.0.0) SFAC = 1.0 IF (RMS.LE.0.0) THEN write(6,*) 'warning: RMS of profile = 0, polynomial fiting', * 'routine will fail' write(6,*) 'have you baselined the profile using base.f ?' ENDIF C C fetching the velocity array and RMS C CALL Recall(Y, Isize, 10 ) C C OK, now go to measure. C IL = PR2 IR = PR3 C exit if invalid measuring method specified; for method 4&5 if C the apropriate fraction is lt 0, check if interactive output C is desired. DO 21 I=1,IMAX PRTLEV = 0 IF (K(I).EQ.0) THEN call putln(" Problem with PR1. Check and try again") GO TO 99 END IF IF((K(I).EQ.4).OR.(K(I).EQ.5))THEN IF(F(I).LT.0.0)THEN 9 WRITE(6,"('lots of output for polynomial methods?')") WRITE(6,"(' y/n [n] or i for interactive')") READ(5,FMT="(A)",ERR=9)RPLY IF(RPLY.EQ.'i')THEN PRTLEV=2 ELSE IF(RPLY.EQ.'y')PRTLEV=1 ENDIF F(I) = -F(I) ENDIF ENDIF IF (K(I).EQ.1) * CALL AREA (X,Y,IL,IR,F(I),SDV(I),JL(I),JR(I),W(I), * V(I),RMS,Pl) IF (K(I).EQ.2) * CALL AREAP (X,Y,IL,IR,F(I),SDV(I),JL(I),JR(I),W(I), * V(I),RMS,Pl) IF(K(I).EQ.3) * CALL AREA2P(X,Y,IL,IR,F(I),SDV(I),JL(I),JR(I),W(I), * V(I),RMS,Pl) IF (K(I).EQ.4) * CALL areaPE(X,Y,Z,T,ISIZE,IL,IR,F(I),SDV(I),JL(I),JR(I), * W(I), V(I),RMS,VRMS,PRTLEV,Pl) IF (K(I).EQ.5) * CALL areaPF(X,Y,Z,T,ISIZE,IL,IR,F(I),SDV(I),JL(I),JR(I), * W(I),V(I),RMS,VRMS,PRTLEV,Pl) WRITE (6,109) IL,IR,SDV(I),V(I),JL(I),JR(I),W(I) 109 FORMAT (" area(",I3,",",I4,") =",E10.4, * " Center =",F11.5,3X," Width(",I3,",",I4,") =",F11.5) C C if the polfit routine has been used then C write the uncertainity in velocity IF((K(I).EQ.4).OR.(K(I).EQ.5))THEN IF(PRTLEV.GE.1)THEN WRITE(6,"('The rms error on V is',F6.1,' km/s')")VRMS ENDIF IF (K(I).EQ.5) THEN VRMS = VRMS*SFAC IF(PRTLEV.GE.1)THEN WRITE(6,"('Scaling by smoothing factor',F5.2,)")SFAC ENDIF ENDIF WRITE(6,"(28X,'DeltaV =',F6.2)")VRMS ENDIF IF (K(I).EQ.1) WD = "M " IF (K(I).EQ.2) WD = "P " IF (K(I).EQ.3) WD = "2P" IF (K(I).EQ.4) WD = "PE" IF (K(I).EQ.5) WD = "PF" WRITE (WORD(I),"(F4.3,A2)") F(I),WD 21 CONTINUE C C C Now write results in STO/RCL 9 C FINT = SDV(1) string = ' area = ' call catd(string, fint, '(E10.4)') call cats(string, ' Ok? (y/n or e (for exit)) ') i = lastc(string) call prompt(string(1:i), 'EXIT_ON_EOF', ans) ans = capitalize(ans) IF (ANS.EQ.'E') GO TO 99 IF (ANS.EQ."N") THEN WRITE (6,'(" Enter desired value of area")') READ (5,*,iostat=iostatus) FINT if (iostatus.ne.0) then call perror(' msr> ') return endif END IF C CALL Recall( Y, Isize, 9 ) Y(139) = SMAX I19 = 140 Y(I19) = FINT C Do 20 I=1,IMAX Y(I19+1) = V(I) Y(I19+2) = W(I) Y(I19+3) = WRD(I) I19 = I19 + 3 20 Continue C Store in STO/RCL 9 the uncertainty in the velocity measurement C VRMS. Meaningful only for the polynomial fitting C velocity subroutine AREAPF or AREAPE; C Store channel numbers at level FR in 241-250, C and slopes of fits to profile sides (method 5) in 251,252 Y(156) = VRMS slope_right = 0. slope_left = 0. do 40, i=1,5 Y(241+(i-1)*2) = JL(i) Y(242+(i-1)*2) = JR(i) if (K(i).eq.5) then slope_left = Z(JL(i)) - Z(JL(i)-1) slope_right= T(JR(i)) - T(JR(i)+1) end if 40 continue Y(251) = slope_left Y(252) = slope_right C CALL Store( Y, Isize, 9 ) C 99 call putln('MSR Done...') RETURN END ******************************************************************************** SUBROUTINE NULLS ( X, PR2, PR3) ******************************************************************************** C Double Precision X(1), Pr2, Pr3 Integer*4 IL, IR, ILL, IRR, JL, JR C IL = PR2 IR = PR3 ILL = IL - 15 IRR = IR + 15 C DO 11 JL = IL,ILL,-1 IF (X(JL).LT.0.) GO TO 20 11 CONTINUE C 20 IL = JL + 1 DO 21 JR = IR,IRR IF (X(JR).LT.0.) GO TO 30 21 CONTINUE C 30 IR = JR - 1 C PR2 = IL PR3 = IR c RETURN END ******************************************************************************** Double Precision FUNCTION SMx ( X, X1, X2 ) ******************************************************************************** C C Finds the maximum between the feature boundaries X1,X2 C Double Precision X(1), x1, x2 Integer*4 J1, J2, J C J1 = X1 J2 = X2 SMX = X(J1) C Do 10 J=J1,J2 SMx = MAX (X(J),SMx) 10 Continue C RETURN END ******************************************************************************** SUBROUTINE AREA (X,Y,IFIRST, LAST, F, SDV, JL, JR, W, V,RMS,Pl) ******************************************************************************** C Double Precision X(1),Y(1),SDV,AVE,BL,FRAC,BR,V,W,FLevel Double Precision Sdv, F,RMS,VR,VL,Pl(1) Integer*4 Ifirst, last, JR, J, JL, I C C C compute a rough estiamte of the area, ignore the C uneven channel spacing in velocity C SDV = 0.0 C Do 10 I=IFIRST,LAST SDV = SDV + X(I) 10 Continue AVE = SDV/dble(LAST-IFIRST+1) C IF(F.GT.1.) THEN WRITE(6,198) 198 FORMAT(" Fraction > 1; do not necessarily " * "expect the correct answer") END IF C FLEVEL = F*AVE C Do 20 I=IFIRST,LAST IF(X(I).GE.FLEVEL) GO TO 55 20 Continue C 55 JL = I BL = FRAC(X(JL-1),X(JL),FLEVEL) VL = Y(JL-1)+BL*(Y(JL)-Y(JL-1)) c Do 30 J = LAST,IFIRST,-1 IF(X(J).GE.FLEVEL) GO TO 66 30 Continue C 66 JR = J BR = FRAC(X(JR),X(JR+1),FLEVEL) VR = Y(JR)+BR*(Y(JR+1)-Y(JR)) C W = VR-VL V = (VR+VL)/2.0 C SDV = 0.0 C C compute area accurately, i.e. using the C fact that the chaanels are not evenly C spaced in velocity C Do 40 I=IFIRST,LAST SDV = SDV + X(I)*(Y(I)-Y(I-1)) 40 Continue C jnc put the channels at which the measurements were made into pl Pl(11) = dble(JL)+BL Pl(12) = dble(JR)+BR 99 RETURN END ******************************************************************************** SUBROUTINE AREAP (X,Y,IFIRST, LAST, F, SDV, JL, JR, W, V,RMS,Pl) ******************************************************************************** C Double Precision X(1),Y(1),SDV, Peak, W, BR,Pl(1), * BL, V, FLevel, F, Frac,RMS, * VL,VR Integer*4 Ifirst, last, Jl, JR , J, I C SDV = 0.0 PEAK = X(IFIRST) C C compute area and peak Do 10 I=IFIRST,LAST IF(X(I).GT.PEAK) PEAK = X(I) SDV = SDV + X(I)*(Y(I)-Y(I-1)) 10 Continue PEAK = PEAK - RMS C IF(F.GT.1.) THEN call putln(" Fraction > 1; do not necessarily "// * "expect the correct answer") END IF C FLEVEL = F*PEAK C Do 20 I=IFIRST,LAST IF(X(I).GE.FLEVEL) GO TO 55 20 Continue C 55 JL = I BL = FRAC(X(JL-1),X(JL),FLEVEL) VL = Y(JL-1)+BL*(Y(JL)-Y(JL-1)) c Do 30 J = LAST,IFIRST,-1 IF(X(J).GE.FLEVEL) GO TO 66 30 Continue c 66 JR = J BR = FRAC(X(JR),X(JR+1),FLEVEL) VR = Y(JR) + BL*(Y(JR+1)-Y(JR)) C W = VR - VL V = (VR +VL)/2.0 C jnc put the channels at which the measurements were made into pl Pl(13) = dble(JL)+BL Pl(14) = dble(JR)+BR RETURN END ******************************************************************************** SUBROUTINE AREA2P (X,Y,IFIRST,LAST,F,SDV,JL,JR,W,V,RMS,Pl) ******************************************************************************** Double Precision X(1),Y(1),Comp,F,Sdv,W,V,FPeak1,Frac,Bl Double Precision FPeak2, Br,RMS,VL,VR,Pl(1) Integer*4 IFirst, Last, k1, k2, Il, I, Jl, JR, IR Integer*4 IFlag, IFF, JJ1, ILL, JJ2 C C Search for correct channel of two horns (IL,IR) C COMP = 0. K1 = IFIRST -10 K2 = IFIRST + 10 C do 10 I=K1,K2 IF(X(I).GE.COMP) THEN COMP = X(I) IL = I END IF 10 Continue COMP = 0. K1 = LAST - 10 K2 = LAST + 10 C Do 20 I=K1,K2 IF(X(I).GE.COMP) THEN COMP = X(I) IR = I END IF 20 Continue C C Find channels at f level of horns and at 1st null C (J1,J2) (JJ1,JJ2) C IFLAG = 1 IFF = IL - 100 C subtract the rms from the peak value FPEAK1 = F*(X(IL)-RMS) C Do 30 I=IL,IFF,-1 IF(I.LT.2) GO TO 201 IF(X(I).LE.FPEAK1.AND.IFLAG.LT.2) THEN JL = I IFLAG = 2 END IF IF(X(I).LE.0.) THEN JJ1 = I + 1 GO TO 201 END IF 30 Continue C 201 IFLAG = 1 ILL = IR + 100 C subtract the rms from the peak value FPEAK2 = F*(X(IR)-RMS) C Do 40 I=IR,ILL IF(I.GT.1005) GO TO 301 IF(X(I).LE.FPEAK2.AND.IFLAG.LT.2) THEN JR = I IFLAG = 2 END IF IF(X(I).LE.0.) THEN JJ2 = I - 1 GO TO 301 END IF 40 Continue C C COMPUTE AREA AND WIDTH 301 SDV = 0.0 Do 50 I=JJ1,JJ2 SDV = SDV + X(I)*(Y(I)-Y(I-1)) 50 Continue C BL = FRAC(X(JL),X(JL+1),FPEAK1) VL = Y(JL) + BL*(Y(JL+1)-Y(JL)) BR = FRAC(X(JR-1),X(JR),FPEAK2) VR = Y(JR-1) + BR*(Y(JR)-Y(JR-1)) W = VR -VL V = (VR + VL)/2.0 C C jnc put the channels at which the measurements were made into pl Pl(15) = dble(JL)+BL Pl(16) = dble(JR)+BR RETURN END *************************************************************************** SUBROUTINE areaPE(X,Y,Z,T,ISIZE,IFIRST,ILAST,F,SDV,IFRL, * IFRR,W,V,RMS,VRMS,PRTLEV,Pl) *************************************************************************** C this subroutine calculates the width and average C velocity of a galaxy spectrum. The two halves of the double C horned profile are treated separately. For each half, the C peak is determined (it is assumed that the two peaks are on either C side of the center of the user specified limits IFIRST and LAST). C Then each on each side of the profile a polynomial (by default C 1st order) is fit to a user specified portion of the profile. C If the fit is poor, (i.e. (chisqr - deg.of.freedom) > 3 times C sqrt(2*deg.of.freedom)), the user is given the option of fitting C a 2nd degree polynomial and/or changing the levels between C which to fit. The FR points (where FR is a user specified fraction C of the peak flux) on either side of the profile are located by C interpolation between the two data points that bracket the C FR point. The interpolaton is done by channel number and then C the velocity corresponding to the interpolated channel number C is computed. The average velocity is taken to be the mean of C these two velocities and the width is taken to be the difference C between these two velocities. C The uncertainity in the velocity measurment is then calculated C from the rms in the baseline and the coefficients of the C fitted polynomial. Finally the area bewteen the channels C IFIRST and ILAST is computed as summation (X(I)*dV(I)), where C dV(I) is the velocity difference between channels I and I-1. C PEAK channels are flagged and the fitted polynomials are plotted. C C the following subroutines are required: C C pksearch: locates the profile horns C plfit : sets up input functions for lfit C lfit : sets up the normal equations matrix,computes the C coefficients of the best fit polynomial C gaussj : does the matrix inversion C covsrt : calculates the covariance matrix of the coefficients C funcs : defines the fitting functions (i.e. polynomials) C INPUT VARIABLES C C X is the array containg the spectrum which is to be measured C Y is the array containing the velocity information for X C ISIZE the length of the arrays C IFIRST is the (user set) channel number of the left edge of the spectrum C ILAST is the (user set) channel number of the right edge of the spectrum C RMS rms noise in the spectrum (as estimated by base.f) C F specifies the levels between which to do the polynomial C fitting and the level at which to measure. A double prec C number of the form 0.abc, where 0.a and 0.b are the upper C and lower fractions of the peak flux between which to C do the fitting and 0.c is the fraction of the peak flux C at which to measure C PRTLEV level of printing to do C OUTPUT VARIABLES C Z fitted polynomial to the left side C T fitted polynomial to the right side C SDV is the area under the spectrum within the IFIRST and ILAST C W is the width of the profile C V is the average velocity of the profile C VRMS is the rms error associated with the velocity measurement C INTERNAL VARIABLES C A(I) parameters of the POLYNOMIAL fit x = A(1) + A(2)*Y + .. C CHAN USED TO CONVERT CHANNEL NUMBERS TO DOUBLE FOR FLAGGING. C COVAR MATRIX CONTAINING COVARIANCE OF THE COEFFICIENTS C OF THE FITTED POLYNOMIALS C dV VELOCITY INTERVAL BETWEEN TWO NEIGHBOURING CHANNELS C FPL LEFT HALF PEAK FLUX LEVEL IN THE SPECTRUM C FPR RIGHT HALF PEAK FLUX LEVEL IN THE SPECTRUM C FR1 FRACTION OF PEAK FLUX AT WHICH TO STOP LINEAR FIT TO PROFILE C FR2 FRACTION OF PEAK FLUX AT WHICH TO START LINEAR FIT TO PROFILE C FR FRACTION OF PEAK FLUX AT WHICH TO CALCULATE WIDTH C I DUMMY INTEGERS FOR LOOPING CONTROL C IFR1L LEFT FR1 CHANNEL C IFR2L LEFT FR2 CHANNEL C IFR1R RIGHT FR1 CHANNEL C IFR2R RIGHT FR1 CHANNEL C IFRL LEFT CHANNEL at level FR C IFRR RIGHT CHANNEL at level FR C IOS STATUS FLAG FOR IO C JUNK VARIABLE IN WHICH TO TRAP WRONG TERMINAL INPUT C RESET FLAGS WETHER USER HAS RESET PARAMETERS C TEMP TEMPORARY VARIABLE USED TO DECODE F C VRMSL RMS ERROR ASSOCIATED WITH THE VELOCITY ESTIMATE C AT FLUX LEVEL FR OF THE LEFT HALF OF THE PROFILE C VRMSR RMS ERROR ASSOCIATED WITH THE VELOCITY ESTIMATE C AT FLUX LEVEL FR OF THE RIGHT HALF OF THE PROFILE C YFRL LEFT FR POINT C YFRR RIGHT FR POINT C DOUBLE PRECISION X(1),Y(1),Z(1),T(1),SDV,V,W,Pl(1) DOUBLE PRECISION FTST1,FTST2,RMS,F Double precision Yvalue, Chan INTEGER*4 IFIRST,ILAST,NORDER,ISIZE, * IFR1L,IFR2L,IFRL,IFR1R,IFR2R, * IFRR,I,IOS,IPL,IPR, * IFRR1,IFRR2,IFRL1,IFRL2,PRTLEV DOUBLE PRECISION VFRL,VFRR,YFRL,YFRR,dV,A(3), * FR1,FR2,FR,FPL,FPR,FPL, * VRMSL,VRMSR,VRMS,COVAR(3,3), * YFRL1,YFRL2,YFRR1,YFRR2,TEMP Double Precision DeltaVF LOGICAL RETRY,ASK CHARACTER JUNK, Key*1 C DECODE F INTO THE APPROPRIATE FRACTIONS c ................................................................... RETRY = .FALSE. 90 IF((.NOT.RETRY).AND.((F.GE.1.0).OR. * F.LE.0.0)) THEN IOS = 0 20 WRITE (6,21) 21 FORMAT (/,'Specify profile levels for fit',/, * 'Check help file for format') READ(5,*,IOSTAT=IOS) F IF(IOS.NE.0) THEN IOS = 0 READ(5,'(a)') JUNK WRITE(6,"('try again...')") GOTO 20 ENDIF IF ((F.GE.1.0).OR.(F.LE.0.0)) THEN WRITE(6,"('try again...')") GOTO 20 ENDIF ENDIF IF (F.NE.0.0) THEN TEMP = F FR1 = DBLE(INT(TEMP*10.0)/10.0) TEMP = TEMP -FR1 FR2 = DBLE(INT(TEMP*100.0)/10.0) TEMP = TEMP - FR2*0.1 FR = DBLE(DNINT(TEMP*1000.0)/10.0) IF((FR.GT.FR1).OR.(FR.LT.FR2).OR. * (FR1.EQ.0.0).OR.(FR2.EQ.0.0).OR. * (FR1.LE.FR2)) THEN WRITE(6,"('Error in polynomial fitting: ')") WRITE(6,"('You need to reset the fitting levels')") GO TO 20 ENDIF ELSE F = FR1 + 0.1*FR2 +0.01*FR ENDIF C LOCATE THE PROFILE HORNS AND FLAG THEM c ..................................................................... CALL PKSEARCH(X,IFIRST,ILAST,IPL,IPR,FPL,FPR) IF(PRTLEV.GE.1)THEN WRITE(6,"('Left peak at channel',i5,' flux',F6.1)")IPL,FPL IF(PRTLEV.GE.2)THEN 151 WRITE(6,"('OK? y/n [y]')") READ(5,FMT="(A)",ERR=151)JUNK IF(JUNK.EQ.'n'.or.JUNK.eq.'N') THEN 152 WRITE(6,"('Flag desired peak channel')") Call Graphmode(.true.) Call Readcursor(key,chan,yvalue) Call Flagpoint(chan) IPL = chan Call Graphmode(.false.) ENDIF ENDIF ENDIF IF(PRTLEV.GE.1) THEN WRITE(6,"('Right peak at channel',i5,' flux',F6.1)")IPR,FPR IF(PRTLEV.GE.2) THEN 153 WRITE(6,"('OK? y/n [y]')") READ(5,FMT="(A)",ERR=153)JUNK IF(JUNK.EQ.'n'.or.JUNK.eq.'N') Then 154 WRITE(6,"('Flag desired peak channel')") Call Graphmode(.true.) Call Readcursor(key,chan,yvalue) Call Flagpoint(chan) IPR = chan Call Graphmode(.false.) ENDIF ENDIF ENDIF C locate the FR1 and FR2 points c ................................................................... DO 102 I=IPL,IFIRST,-1 IF(X(I).LT.(FR1*FPL)) GOTO 201 102 CONTINUE 201 IFR1L = I DO 103 I= IPL,IFIRST,-1 IF(X(I).LT.FR2*FPL) GOTO 202 103 CONTINUE 202 IFR2L = I IF(PRTLEV.GE.1)THEN WRITE(6,"('We fit left side between chan',i5,i5)")IFR1L,IFR2L IF(PRTLEV.GE.2)THEN 155 WRITE(6,"('OK? y/n [y]')") READ(5,FMT="(A)",ERR=155)JUNK IF(JUNK.EQ.'n'.or.JUNK.eq.'N')THEN 156 WRITE(6,"('Flag the desired 1st fit channel')") Call Graphmode(.true.) Call Readcursor(key,chan,yvalue) Call Flagpoint(chan) IFR1L = chan Call Graphmode(.false.) WRITE(6,"('And now flag the desired last fit channel')") Call Graphmode(.true.) Call Readcursor(key,chan,yvalue) Call Flagpoint(chan) IFR2L = chan Call Graphmode(.false.) ENDIF ENDIF ENDIF IF(IFR1L.LE.IFR2L) THEN IFR1L = IFR2L + 1 WRITE (6,"('**** Warning: upper and lower profile ', * 'points coincide!')") ENDIF DO 104 I=IPR, ILAST IF(X(I).LT.FR1*FPR) GOTO 203 104 CONTINUE 203 IFR1R = I DO 105 I=IPR,ILAST IF(X(I).LT.FR2*FPR) GOTO 204 105 CONTINUE 204 IFR2R = I IF(PRTLEV.GE.1)THEN WRITE(6,"('Fit right side between chan',I5,I5)")IFR1R,IFR2R IF(PRTLEV.GE.2)THEN 157 WRITE(6,"('OK? y/n [y]')") READ(5,FMT="(A)",ERR=157)JUNK IF(JUNK.EQ.'n'.or.junk.eq.'N')THEN 158 WRITE(6,"('Flag desired 1st fit channel')") Call Graphmode(.true.) Call Readcursor(key,chan,yvalue) Call Flagpoint(chan) IFR1R = chan Call Graphmode(.false.) WRITE(6,"('And now flag the desired last fit channel')") Call Graphmode(.true.) Call Readcursor(key,chan,yvalue) Call Flagpoint(chan) IFR2R = chan Call Graphmode(.false.) ENDIF ENDIF ENDIF IF(IFR1R.GE.IFR2R) THEN IFR2R = IFR1R + 1 WRITE (6,"('**** Warning: upper and lower profile ', * 'points coincide!')") ENDIF IF(.NOT.RETRY) NORDER = 1 C reduce the peak flux levels by the RMS C FPL = FPL - RMS C FPR = FPR - RMS C locate the FR points c .............................................................. IFRL1 = IFIRST DO 601 I = IPL, IFIRST,-1 IF(X(I).LT.FR*FPL) GOTO 602 601 CONTINUE 602 IFRL1 = I CALL PLFIT(X,IFRL1,IFRL1+1,RMS,1,A,FTST1,COVAR) YFRL1 = (FPL*FR -A(1))/A(2) IFRL2 = IPL DO 603 I = IFIRST,IPL IF(X(I).GT.FR*FPL) GOTO 604 603 CONTINUE 604 IFRL2 = I CALL PLFIT(X,IFRL2-1,IFRL2,RMS,1,A,FTST1,COVAR) YFRL2 = (FPL*FR -A(1))/A(2) YFRL = (YFRL1+YFRL2)/2.0 IFRL = INT(YFRL) IFRR1 = IPR DO 605 I = IPR, ILAST IF(X(I).LT.FR*FPR) GOTO 606 605 CONTINUE 606 IFRR1 = I CALL PLFIT (X,IFRR1-1,IFRR1,RMS,1,A,FTST2,COVAR) YFRR1 = (FPR*FR -A(1))/A(2) IFRR2 = IPR DO 607 I = ILAST,IPR,-1 IF(X(I).GT.FR*FPR) GOTO 608 607 CONTINUE 608 IFRR2 = I CALL PLFIT (X,IFRR2,IFRR2+1,RMS,1,A,FTST2,COVAR) YFRR2 = (FPR*FR -A(1))/A(2) YFRR = (YFRR1 + YFRR2)/2.0 IFRR =INT(YFRR) C fit a polynomial of order NORDER to the sides C of the profile and calculate the velocity C uncertainity and c .............................................................. CALL PLFIT (X,IFR2L,IFR1L,RMS,NORDER,A,FTST1,COVAR) IF (NORDER.EQ.1) THEN VRMSL = RMS/A(2)*(Y(IFRL)- Y(IFRL-1)) DO 501 I=1,ISIZE Z(I)=A(1)+A(2)*I 501 CONTINUE ENDIF IF(NORDER.EQ.2) THEN VRMSL = RMS/(2.0*A(3)*YFRL+A(2))*(Y(IFRL)-Y(IFRL-1)) DO 502 I=1,ISIZE Z(I)=A(1)+A(2)*I+A(3)*I*I 502 CONTINUE ENDIF c prevent the velocity error on left side from exceeding c the velocity range of the whole fit. DeltaVF = Abs(Y(IFR2L)-Y(IFR1L)) C C TYPO!!!!!! VRSML ---->>>> VRmsL If(VRmsL.gt.DeltaVF) VRMSL = DeltaVF IF(PRTLEV.GE.1)THEN WRITE(6,"('Left half rms in channel number units',F7.2)") * VRMSL/(Y(IFRL)-Y(IFRL-1)) ENDIF CALL PLFIT (X,IFR1R,IFR2R,RMS,NORDER,A,FTST2,COVAR) IF (NORDER.EQ.1) THEN VRMSR = RMS/A(2)*(Y(IFRR)-Y(IFRR-1)) DO 503 I=1,ISIZE T(I)=A(1)+A(2)*I 503 CONTINUE ENDIF IF(NORDER.EQ.2) THEN VRMSR = RMS/(2.0*A(3)*YFRR+A(2))*(Y(IFRR)-Y(IFRR-1)) DO 504 I=1,ISIZE T(I)=A(1)+A(2)*I+A(3)*I*I 504 CONTINUE ENDIF c prevent the velocity error on right side from exceeding c the velocity range of the whole fit. DeltaVF = Abs(Y(IFR1R)-Y(IFR2R)) If (VRMSR.gt.DeltaVF) VRMSR = DeltaVF IF(PRTLEV.GE.1)THEN WRITE(6,"('Right half rms in channel number units',F7.2)") * VRMSR/(Y(IFRR)-Y(IFRR-1)) ENDIF C if the fit is poor let the user reset the parameters c ................................................................ IF ((DABS(FTST1).GT.3.0).OR.(DABS(FTST2).GT.3.0)) THEN WRITE(6,10) FTST1,FTST2,IFR1L-IFR2L+1,IFR2R-IFR1R+1 10 FORMAT('Poor fit to profile; the normalized ch1sqr are', * /,f6.1,f61, /, * '# of points used to fit l & r halves are', * x,i3,x,i3 ) RETRY = ASK('Fit a 2nd degree polynomial or change levels? ' * , ' y/n [y]' ,.TRUE.) IF (RETRY) THEN 22 WRITE(6,23) 23 FORMAT(/,'Specify polynomial order (le.2) and profile levels', * 'between which to fit, in standard input format',/, * '[e.g. 2 0.825; ... or for default enter', * '0 0 ==> use 2nd order & original profile levels]') READ(5,*,IOSTAT=IOS) NORDER,F IF(IOS.NE.0) THEN IOS = 0 READ (5,'(a)') JUNK WRITE(6,"( 'try again...')") GOTO 22 ENDIF IF ((F.GE.1.0).OR.(F.LT.0.0).OR. * (NORDER.LT.0).OR.(NORDER.GT.2)) THEN WRITE(6,"('try again...')") GOTO 22 ENDIF IF(NORDER.LE.0) NORDER = 2 GOTO 90 ENDIF ENDIF C calculate the velocity width and average velocity C .................................................................. IFRL = INT(YFRL) VFRL = Y(IFRL) + (YFRL -DBLE(IFRL))*(Y(IFRL+1)-Y(IFRL)) IFRR = INT(YFRR) VFRR = Y(IFRR) + (YFRR -DBLE(IFRR))*(Y(IFRR+1)-Y(IFRR)) V = (VFRL + VFRR)/2.0 W = (VFRR - VFRL) VRMS = SQRT((VRMSL**2 +VRMSR**2)/2.0) C calculate the integrated flux in bounded region C ................................................................... SDV = 0.0 DO 106 I = IFIRST, ILAST dV = Y(I) - Y(I-1) SDV = SDV +X(I)*dV 106 CONTINUE C return the peak and fit channel numbers in the pl array C .................................................................... pl(17) = dble(IPL) pl(18) = dble(IPR) pl(19) = YFRL pl(20) = YFRR pl(31) = dble(IFR1L) pl(32) = dble(IFR2L) pl(33) = dble(IFR1R) pl(34) = dble(IFR2R) RETURN END *************************************************************************** SUBROUTINE areaPF(X,Y,Z,T,ISIZE,IFIRST,ILAST,F,SDV, * IFRL,IFRR,W,V,RMS,VRMS,PRTLEV,Pl) *************************************************************************** C this subroutine calculates the width and average C velocity of a galaxy spectrum. The two halves of the double C horned profile are treated separately. For each half, the C peak is determined (it is assumed that the two peaks are on either C side of the center of the user specified limits IFIRST and LAST). C Then each on each side of the profile a polynomial (by default C 1st order) is fit to a user specified portion of the profile. C If the fit is poor, (i.e. (chisqr - deg.of.freedom) > 3 times C sqrt(2*deg.of.freedom)), the user is given the option of fitting C a 2nd degree polynomial and/or changing the levels between C which to fit. Finaly using the fitted polynomial, the FR points C (where FR is a user specified fraction of the peak flux) on C either side of the profile is located. The interpolaton C is done by channel number and then the velocity corresponding C to the interpolated channel number is computed. The average velocity C is taken to be the mean of these two velocities and the C width is taken to be the difference bwteen these two velocities. C The uncertainity in the velocity measurment is then calculated C from the covariance matrix of the coefficients of the fitted C polynomial. Finally the area of the profile between the limits C IFIRST and ILAST is computed as summation (X(I)*dV(I)), where C dV(I) is the velocity difference between channels I and I-1. C PEAK channels are flagged and the fitted polynomials are plotted. C the following subroutines are required: C pksearch : locates the profile horns C plfit : sets up input functions for lfit C lfit : sets up the normal equations matrix,computes the C coefficients of the best fit polynomial C gaussj : does the matrix inversion C covsrt : calculates the covariance matrix of the coefficients C funcs : defines the fitting functions (i.e. polynomials) C C INPUT VARIABLES C X is the array containg the spectrum which is to be measured C Y is the array containing the velocity information for X C ISIZE the length of the array C IFIRST is the (user set) channel number of the left edge of the spectrum C ILAST is the (user set) channel number of the right edge of the spectrum C RMS rms noise in the spectrum (as estimated by base.f) C F specifies the levels between which to do the polynomial C fitting and the level at which to measure. A double prec C number of the form 0.abc, where 0.a and 0.b are the upper C and lower fractions of the peak flux between which to C do the fitting and 0.c is the fraction of the peak flux C at which to measure C PRTLEV level of printing to do C OUTPUT VARIABLES C Z fitted polynomial to the left hand side C T fitted polynomial to the right hand side C SDV is the area under the spectrum within the IFIRST and ILAST C W is the width of the profile C V is the average velocity of the profile C VRMS is the rms error associated with the velocity measurement C INTERNAL VARIABLES C A(I) parameters of the POLYNOMIAL fit x = A(1) + A(2)*Y + .. C B(I) DUMMY VECTOR USED TO CONDENSE RMS CALCULATION CODE C CHAN USED TO CONVERT CHANNEL NUMBERS TO DOUBLE FOR FLAGGING. C COVAR MATRIX CONTAINING COVARIANCE OF THE COEFFICIENTS C OF THE FITTED POLYNOMIALS C dV VELOCITY INTERVAL BETWEEN TWO NEIGHBOURING CHANNELS C FPL LEFT HALF PEAK FLUX LEVEL IN THE SPECTRUM C FPR RIGHT HALF PEAK FLUX LEVEL IN THE SPECTRUM C FR1 FRACTION OF PEAK FLUX AT WHICH TO STOP LINEAR FIT TO PROFILE C FR2 FRACTION OF PEAK FLUX AT WHICH TO START LINEAR FIT TO PROFILE C FR FRACTION OF PEAK FLUX AT WHICH TO CALCULATE WIDTH C FRMSL RMS ERROR OF THE FLUX AT THE INTERPOLATED LEVEL FR OF C THE PEAK FLUX C I,J DUMMY INTEGERS FOR LOOPING CONTROL C IFR1L LEFT FR1 CHANNEL C IFR2L LEFT FR2 CHANNEL C IFR1R RIGHT FR1 CHANNEL C IFR2R RIGHT FR1 CHANNEL C IFRL LEFT CHANNEL at FR level C IFRR RIGHT CHANNEL at FR level C IOS STATUS FLAG FOR IO C JUNK VARIABLE IN WHICH TO TRAP WRONG TERMINAL INPUT C RESET FLAGS WETHER USER HAS RESET PARAMETERS C TEMP TEMPORARY VARIABLE USED TO DECODE F C VRMSL RMS ERROR ASSOCIATED WITH THE VELOCITY ESTIMATE C AT FLUX LEVEL FR OF THE LEFT HALF OF THE PROFILE C VRMSR RMS ERROR ASSOCIATED WITH THE VELOCITY ESTIMATE C AT FLUX LEVEL FR OF THE RIGHT HALF OF THE PROFILE C YFRL LEFT FR POINT C YFRR RIGHT FR POINT DOUBLE PRECISION X(1),Y(1),Z(1),T(1),Pl(1),SDV,V,W Double precision Yvalue, Chan DOUBLE PRECISION FTST1,FTST2,RMS ,F INTEGER*4 IFIRST,ILAST,NORDER,ISIZE, * IFR1L,IFR2L,IFRL,IFR1R,IFR2R, * IFRR,I,IOS,IPL,IPR,J,PRTLEV DOUBLE PRECISION VFRL,VFRR,YFRL,YFRR,dV,A(3), * FR1,FR2,FR,FPL,FPR,FPL,B(3), * FRMS,VRMSL,VRMSR,VRMS, * COVAR(3,3),TEMP Double Precision Epsl_ch,Epsr_ch,DR_max,DL_max LOGICAL RETRY,ASK CHARACTER JUNK, Key*1 C DECODE F INTO THE APPROPRIATE FRACTIONS C .................................................................... RETRY = .FALSE. 90 IF((.NOT.RETRY).AND.((F.GE.1.0).OR. * F.LE.0.0)) THEN IOS = 0 20 WRITE (6,21) 21 FORMAT (/,'Specify profile levels for fit',/, * 'check help file for format') READ(5,*,IOSTAT=IOS) F IF(IOS.NE.0) THEN IOS = 0 READ(5,'(a)') JUNK WRITE(6,"('try again...')") GOTO 20 ENDIF IF ((F.GE.1.0).OR.(F.LE.0.0)) THEN WRITE(6,"( 'try again...')") GOTO 20 ENDIF ENDIF IF (F.NE.0.0) THEN TEMP = F FR1 = DBLE(INT(TEMP*10.0)/10.0) TEMP = TEMP -FR1 FR2 = DBLE(INT(TEMP*100.0)/10.0) TEMP = TEMP - FR2*0.1 FR = DBLE(DNINT(TEMP*1000.0)/10.0) IF((FR.GT.FR1).OR.(FR.LT.FR2).OR. * (FR1.EQ.0.0).OR.(FR2.EQ.0.0).OR. * (FR1.LE.FR2)) THEN WRITE(6,"('Error in polynomial fitting: ')") WRITE(6,"('You need to reset the fitting levels')") GO TO 20 ENDIF ELSE F = FR1+0.1*FR2+0.01*FR ENDIF C LOCATE THE PROFILE HORNS AND FLAG THEM C ................................................................. CALL PKSEARCH(X,IFIRST,ILAST,IPL,IPR,FPL,FPR) IF(PRTLEV.GE.1)THEN WRITE(6,"('Left peak at channel',i5,' flux',f6.1)")IPL,FPL IF(PRTLEV.GE.2)THEN 151 WRITE(6,"('OK? y/n [y]')") READ(5,FMT="(A)",ERR=151)JUNK IF(JUNK.EQ.'n'.or.JUNK.eq.'N')THEN 152 WRITE(6,"('Flag desired peak channel')") Call Graphmode(.true.) Call Readcursor(key,chan,yvalue) Call Flagpoint(chan) IPL = chan Call Graphmode(.false.) ENDIF ENDIF ENDIF IF(PRTLEV.GE.1)THEN WRITE(6,"('Right peak at channel',i5,' flux',F6.1)")IPR,FPR IF(PRTLEV.GE.2)THEN 153 WRITE(6,"('OK? y/n [y]')") READ(5,FMT="(A)",ERR=153)JUNK IF(JUNK.EQ.'n'.or.JUNK.eq.'N')THEN 154 WRITE(6,"('Flag desired peak channel')") Call Graphmode(.true.) Call Readcursor(key,chan,yvalue) Call Flagpoint(chan) IPR = chan Call Graphmode(.false.) READ(5,"(I)",ERR=154)IPR ENDIF ENDIF ENDIF C locate the FR1 and FR2 points C ............................................................. DO 102 I=IPL,IFIRST,-1 IF(X(I).LT.(FR1*FPL)) GOTO 201 102 CONTINUE 201 IFR1L = I DO 103 I= IPL,IFIRST,-1 IF(X(I).LT.FR2*FPL) GOTO 202 103 CONTINUE 202 IFR2L = I IF(PRTLEV.GE.1)THEN WRITE(6,"('Fit left side between chan',I5,I5)")IFR2L,IFR1L IF(PRTLEV.GE.2)THEN 155 WRITE(6,"('OK? y/n [y]')") READ(5,FMT="(A)",ERR=155)JUNK IF(JUNK.EQ.'n'.or.JUNK.eq.'N')THEN 156 WRITE(6,"('Flag 1st desired fit channel')") Call Graphmode(.true.) Call Readcursor(key,chan,yvalue) Call Flagpoint(chan) IFR2L = chan Call Graphmode(.false.) WRITE(6,"('Now flag last desired fit channel')") Call Graphmode(.true.) Call Readcursor(key,chan,yvalue) Call Flagpoint(chan) IFR1L = chan Call Graphmode(.false.) ENDIF ENDIF ENDIF IF(IFR1L.LE.IFR2L) THEN IFR1L = IFR2L + 1 WRITE (6,"('**** Warning: upper and lower profile ', * 'points coincide!')") ENDIF DO 104 I=IPR, ILAST IF(X(I).LT.FR1*FPR) GOTO 203 104 CONTINUE 203 IFR1R = I DO 105 I= IPR,ILAST IF(X(I).LT.FR2*FPR) GOTO 204 105 CONTINUE 204 IFR2R = I IF(PRTLEV.GE.1)THEN WRITE(6,"('Fit right side between chan',I5,I5)")IFR1R,IFR2R IF(PRTLEV.GE.2)THEN 157 WRITE(6,"('OK? y/n [y]')") READ(5,FMT="(A)",ERR=157)JUNK IF(JUNK.EQ.'n'.or.JUNK.eq.'N')THEN 158 WRITE(6,"('Flag 1st desired fit channel')") Call Graphmode(.true.) Call Readcursor(key,chan,yvalue) Call Flagpoint(chan) IFR1R = chan Call Graphmode(.false.) WRITE(6,"('Now flag last desired fit channel')") Call Graphmode(.true.) Call Readcursor(key,chan,yvalue) Call Flagpoint(chan) IFR2R = chan Call Graphmode(.false.) ENDIF ENDIF ENDIF IF(IFR1R.GE.IFR2R) THEN IFR2R = IFR1R + 1 WRITE (6,"('****Warning: upper and lower profile ', * 'points coincide')") ENDIF IF(.NOT.RETRY) NORDER = 1 C reduce the peak flux levels by the RMS C FPL = FPL - RMS C FPR = FPR - RMS C fit a polynomial of order NORDER to the left half C of the profile and locate the FR point.store fit in C Z array. C ........................................................... DL_MAX = Real(IFR1L) - Real(IFR2L) CALL PLFIT (X,IFR2L,IFR1L,RMS,NORDER,A,FTST1,COVAR) IF (NORDER.EQ.1) THEN YFRL = (FPL*FR -A(1))/A(2) IFRL = INT(YFRL) ENDIF DO 501 I= 1,ISIZE Z(I) = A(1) +A(2)*I 501 CONTINUE IF(NORDER.EQ.2) THEN YFRL = (-A(2) + SQRT(A(2)*A(2) -4.0*(A(1)-FR*FPL)*A(3))) YFRL = YFRL /(2.0*A(3)) IFRL = INT(YFRL) IF((IFRL.LT.IFR2L).OR.(IFRL.GT.IFR1L)) THEN YFRL = (-A(2)- SQRT(A(2)*A(2)-4.0*(A(1)-FR*FPL)*A(3))) YFRL = YFRL /(2.0*A(3)) ENDIF ENDIF DO 502 I= 1,ISIZE Z(I)=A(1)+A(2)*I+A(3)*I*I 502 CONTINUE C Calculate the uncertainity in the velocity measurement C and constrain it to less than the interval of fit C .............................................................. B(1) = 1.0 B(2) = YFRL B(3) = YFRL*YFRL FRMS = 0.0 DO 301 I = 1, NORDER + 1 DO 302 J = 1, NORDER + 1 FRMS = FRMS + B(I)*B(J)*COVAR(I,J) 302 CONTINUE 301 CONTINUE If (NORDER.eq.1) EPSL_CH = FRMS/A(2) If (NORDER.eq.2) EPSL_CH = FRMS/(2.0*A(3)+YFRL+A(2)) If (EPSL_CH.le.DL_MAX) then VRMSL = EPSL_CH*(Y(IFRL)-Y(IFRL-1)) else VRMSL = DL_MAX*(Y(IFRL)-Y(IFRL-1)) End if IF(PRTLEV.GE.1)THEN WRITE(6,"('Left half rms in channel number units',F7.2)") * VRMSL/(Y(IFRL)-Y(IFRL-1)) ENDIF C Fit a polynomial of order NORDER to the right half C of the profile and locate the FR point. Store fit in C T array. C .............................................................. DR_MAX = Real(IFR2R) - Real(IFR1R) CALL PLFIT (X,IFR1R,IFR2R,RMS,NORDER,A,FTST2,COVAR) IF (NORDER.EQ.1) THEN YFRR = (FPR*FR -A(1))/A(2) IFRR = INT(YFRR) ENDIF DO 503 I=1,ISIZE T(I)=A(1)+A(2)*I 503 CONTINUE IF(NORDER.EQ.2) THEN YFRR = (-A(2) + SQRT(A(2)*A(2) -4.0*(A(1)-FR*FPR)*A(3))) YFRR = YFRR /(2.0*A(3)) IFRR = INT(YFRR) IF((IFRR.LT.IFR1R).OR.(IFRR.GT.IFR2R)) THEN YFRR =(-A(2)-SQRT(A(2)*A(2)-4.0*(A(1)-FR*FPR)*A(3))) YFRR = YFRR /(2.0*A(3)) ENDIF ENDIF DO 504 I=1,ISIZE T(I)=A(1)+A(2)*I+A(3)*I*I 504 CONTINUE C Calculate the uncertainity in the velocity measurement C and constrain it to be less than the interval of fit C ................................................................ B(1) = 1.0 B(2) = YFRR B(3) = YFRR*YFRR FRMS = 0.0 DO 303 I = 1, NORDER + 1 DO 304 J = 1, NORDER + 1 FRMS = FRMS + B(I)*B(J)*COVAR(I,J) 304 CONTINUE 303 CONTINUE If (NORDER.eq.1) EPSR_CH = FRMS/A(2) If (NORDER.eq.2) EPSR_CH = FRMS/(2.0*A(3)+YFRR+A(2)) If (EPSR_CH.le.DR_MAX) then VRMSR = EPSR_CH*(Y(IFRR)-Y(IFRR-1)) else VRMSR = DL_MAX*(Y(IFRR)-Y(IFRR-1)) End if IF(PRTLEV.GE.1)THEN WRITE(6,"('Right half rms in channel number units',F7.2)") * VRMSR/(Y(IFRR)-Y(IFRR-1)) ENDIF C If the fit is poor let the user reset the parameters C ................................................................. IF ((DABS(FTST1).GT.3.0).OR.(DABS(FTST2).GT.3.0)) THEN WRITE(6,10) FTST1,FTST2,IFR1L-IFR2L+1,IFR2R-IFR1R+1 10 FORMAT('Poor fit to profile; the normalized ch1sqr are', * /,f6.1,f61, /, * '# of points used to fit l & r halves are', * x,i3,x,i3 ) RETRY = ASK('Fit a 2nd degree polynomial or change levels? ' * , ' y/n [y]' ,.TRUE.) IF (RETRY) THEN 22 WRITE(6,23) 23 FORMAT(/,'Specify polynomial order (le.2) and profile levels', * 'between which to fit, in standard input format',/, * '[e.g. 2 0.825; ... or for default enter', * '0 0 ==> use 2nd order & original profile levels]') READ(5,*,IOSTAT=IOS) NORDER,F IF(IOS.NE.0) THEN IOS = 0 READ (5,'(a)') JUNK WRITE(6,"('try again...')") GOTO 22 ENDIF IF ((F.GE.1.0).OR.(F.LT.0.0).OR. * (NORDER.LT.0).OR.(NORDER.GT.2)) THEN WRITE(6,"('try again...')") GOTO 22 ENDIF IF(NORDER.LE.0) NORDER = 2 GOTO 90 END IF END IF C calculate the velocity width and average velocity C ................................................................. IFRL = INT(YFRL) VFRL = Y(IFRL) + (YFRL -DBLE(IFRL))*(Y(IFRL+1)-Y(IFRL)) IFRR = INT(YFRR) VFRR = Y(IFRR) + (YFRR -DBLE(IFRR))*(Y(IFRR+1)-Y(IFRR)) V = (VFRL + VFRR)/2.0 W = (VFRR - VFRL) VRMS = SQRT((VRMSL**2 +VRMSR**2)/2.0) C calculate the integrated flux in bounded region C .................................................................. SDV = 0.0 DO 106 I = IFIRST, ILAST dV = Y(I) - Y(I-1) SDV = SDV +X(I)*dV 106 CONTINUE C return the peak and fit channel numbers in the pl array C .................................................................. pl(21) = dble(IPL) pl(22) = dble(IPR) pl(23) = YFRL pl(24) = YFRR pl(35) = dble(IFR1L) pl(36) = dble(IFR2L) pl(37) = dble(IFR1R) pl(38) = dble(IFR2R) RETURN END ************************************************************************* subroutine pksearch(X,IFIRST,ILAST,IPL,IPR,FPL,FPR) ************************************************************************* C this subroutine is a first step at automating the velocity C measurement for double horned profiles. The subroutine C will return the fluxes and channel numbers of the two C horns. First the entire profile between the user flaged C limits IFIRST and ILAST (IFIRST < ILAST), is searched C for peaks, a peak being defined as a channel which has C more flux than NISO channels on either side of it. NISO C is scaled linearly with the width of the profile (i.e. C IFIRST - ILAST) and is 5 when (IFIRST - ILAST) = 40). C If the profile has only one peak, then it is treated C as a double horned porfile in which the peaks are C coincident. If there are more than two peaks then C the outer two peaks are regarded as the horns of C profile UNLESS one of the inner peaks is more than C a factor of FACT times larger than an exterior peak. C in that case the closest exterior peak is replaced C by the large interior peak. In case there is more than C one interior peak that is a factor of FACT larger than C the exterior peak, the outer most large peak is C used. C C INPUT VARIABLES C C X array containg the profile C IFIRST user defined left end of the profile C ILAST user defined right end of the profile C C C OUTPUT VARIABLES C C IPL channel number of left peak C IPR channel number of right peak C FPL Flux of the left peak C FPR Flux of the right peak C C INTERNAL VARIABLES C C PMAX max number of peaks allowed C NPEAK number of peaks in the profile C PKCHAN array containing channel numbers of the peaks C NISO a peak is defined as a channel with more C flux than NISO channels on either side C ICNTR the midpoint of the profile C PCNTR the channel number of the central peak C I,J,K dummy integers for looping control C FACT an internal peak has to be greater than C FACT times an exterior peak to be considered C the profile horn C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER*4 PMAX PARAMETER(PMAX=20) DOUBLE PRECISION X(1),FPL,FPR INTEGER*4 IPR,IPL,IFIRST,ILAST INTEGER*4 NISO,NPEAK,PKCHAN(PMAX),PCNTR,ICNTR, * I,J,K DOUBLE PRECISION FACT IPL = IFIRST FPL = X(IFIRST) IPR = ILAST FPR = X(ILAST) FACT = 2.0 C C A CHANNEL HAS TO BE LARGER IN VALUE THAN NISO C CHNNELS ON EITHER SIDE TO BE CONSIDERED A C PEAK C NISO = INT(5.0*(ILAST - IFIRST)/40.0) IF(NISO.LT.3) NISO = 3 if(niso.gt.7) niso = 7 C FIND ALL PEAKS BETWEEN THE FLAGED LIMITS K=0 DO 390 I=IFIRST,IFIRST+NISO DO 391 J=1,NISO IF(X(I).LT.X(I+J)) GOTO 393 391 CONTINUE DO 392 J=1,I-IFIRST+1 IF(X(I).LT.X(I-J+1)) GOTO 393 392 CONTINUE K=K+1 PKCHAN(K)=I 393 CONTINUE 390 CONTINUE DO 400 I =IFIRST+NISO, ILAST-NISO DO 401 J = 1,NISO IF((X(I).LT.X(I+J)).OR.(X(I).LT.X(I-J))) GOTO 402 401 CONTINUE K = K + 1 PKCHAN(K) = I 402 CONTINUE 400 CONTINUE DO 394 I=ILAST-NISO,ILAST DO 395 J=1,NISO IF(X(I).LT.X(I-J)) GOTO 397 395 CONTINUE DO 396 J=1,ILAST-I+1 IF(X(I).LT.X(I+J-1)) GOTO 397 396 CONTINUE K=K+1 PKCHAN(K)=I 397 CONTINUE 394 CONTINUE NPEAK = K C C IF THERE IS ONLY ONE PEAK, THE LEFT AND RIGHT C SIDE HORNS ARE COINCIDENT C IF (NPEAK.EQ.1) THEN IPL = PKCHAN(1) FPL = X(IPL) IPR = IPL FPR = FPL ELSE C C IF THERE ARE ONLY TWO PEAKS, THEN STORE THE C PEAK FLUXES AND CHANNEL NUMBERS C IF(NPEAK.EQ.2) THEN IPL = PKCHAN(1) IPR = PKCHAN(2) FPL = X(IPL) FPR = X(IPR) ELSE C C IF THERE ARE MORE THAN TWO PEAKS TREAT THE C OUTERMOST TWO PEAKS AS THE PROFILE HORN, C UNLESS ONE OF THE INTERIOR PEAKS IS MORE C THAN FACT TIMES LARGER THAN THE EXTERIOR C PEAKS (SEE EXPLANATION AT START) C IPL = PKCHAN(1) FPL = X(IPL) IPR = PKCHAN(NPEAK) FPR = X(IPR) ICNTR = (PKCHAN(NPEAK) -PKCHAN(1))/2 PCNTR = 2 DO 403, J = 2,NPEAK-1 IF(PKCHAN(J).LT.ICNTR) PCNTR = J 403 CONTINUE J = PCNTR 404 IF(J.GE.2) THEN IF(X(PKCHAN(J)).GT.FACT*X(PKCHAN(1))) THEN IPL = PKCHAN(J) FPL = X(IPL) ENDIF J = J -1 GOTO 404 ENDIF J = PCNTR + 1 405 IF(J.LE.NPEAK -1) THEN IF(X(PKCHAN(J)).GT.FACT*X(PKCHAN(NPEAK))) THEN IPR = PKCHAN(J) FPR = X(IPR) ENDIF J = J + 1 GOTO 405 ENDIF ENDIF ENDIF RETURN END ******************************************************************************* SUBROUTINE PLFIT(X,I1,I2,RMS,NORDER,A,FITRMS,COVAR) **************************************************************************** CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C sets up the input parameters for lfit, computes goodness of C fit C C INPUT PARAMETERS C C X array of values for fitting the polynomial C I1,I2, limits between which to fit C RMS rms error in the data points X C NORDER fitting polynomial order C COVAR covaraince matrix of coefficients A(I) C C OUTPUT PARAMETERS C C A vector containing coeficients of the fit polynomial C X = A(1) +A(2)*Y +.. C FITRMS goodness of fit normalised to 0 mean, unit variance C C INTERNAL VARIABLES C C XA X shifted left by I1 C YA array containg the channel numbers from I1 to I2, C this is the independent variable for the fitting C routine C NMAX max number of data points allowed C NCVM max number of fitting coefficients allowed C SIG array containg the rms of each data point C i.e. SIG(I) = RMS for all I C CHISQR chisqr of the fit C MFIT number of coeficients to fit C MA number of coeficients (i.e. MA - MFIT coefficients C are held fixed C NDATA number of data points C NFREE number of degrees of freedom C I dummy integer for looping control C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC integer*4 NMAX,NCVM PARAMETER (NMAX = 50) PARAMETER (NCVM = 3) DOUBLE PRECISION X(1),A(NCVM),FITRMS,RMS INTEGER*4 I1,I2,NORDER DOUBLE PRECISION XA(NMAX),YA(NMAX),SIG(NMAX),A(NCVM) DOUBLE PRECISION COVAR(NCVM,NCVM),CHISQR INTEGER*4 MFIT,MA,NDATA,LISTA(NCVM),I,NFREE NDATA = I2 -I1 + 1 DO 100 I = 1,NDATA XA(I) = X(I+I1 -1) YA(I) = DBLE(I+I1 -1) SIG(I) = RMS 100 CONTINUE MFIT = NORDER + 1 MA = MFIT DO 101 I =1, MA LISTA(I) = 0 IF(I.LE.(NORDER+1)) LISTA(I) = I 101 CONTINUE CALL LFIT(YA,XA,SIG,NDATA,A,MA,LISTA,MFIT,COVAR,NCVM,CHISQR) NFREE = NDATA - MFIT IF (NFREE.GT.0) THEN FITRMS = (CHISQR - NFREE)/SQRT(2.0*NFREE) ELSE FITRMS = 0.0 IF(NFREE.LT.0) THEN WRITE(6,20)NORDER,NDATA 20 FORMAT('warning:',i3,'order polynomial fit to', * i3,'points') ENDIF ENDIF RETURN END ****************************************************************************** SUBROUTINE FUNCS (X,AFUNC,MA) ****************************************************************************** INTEGER*4 MA DOUBLE PRECISION AFUNC(MA),X INTEGER*4 I AFUNC(1) = 1.0 DO 100 I = 2 ,MA AFUNC(I) = AFUNC(I-1)*X 100 CONTINUE RETURN END *************************************************************************** SUBROUTINE LFIT(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,COVAR,NCVM,CHISQ) ************************************************************************** CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C numerical recipes (Press et.al. ) rotuine C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC integer*4 MMAX PARAMETER (MMAX=50) DOUBLE PRECISION X(NDATA),Y(NDATA),SIG(NDATA), * COVAR(NCVM,NCVM),AFUNC(MMAX), * BETA(MMAX),CHISQ,SIG2I,SUM, * A(MA),YM,WT INTEGER*4 LISTA(MA),NDATA,NCVM,MA,MFIT,I,KK, * IHIT,K,J KK=MFIT+1 DO 12 J=1,MA IHIT=0 DO 11 K=1,MFIT IF (LISTA(K).EQ.J) IHIT=IHIT+1 11 CONTINUE IF (IHIT.EQ.0) THEN LISTA(KK)=J KK=KK+1 ELSE IF (IHIT.GT.1) THEN write(6,*) 'Improper set in LISTA' goto 25 ENDIF 12 CONTINUE IF (KK.NE.(MA+1)) then write(6,*) 'Improper set in LISTA' goto 25 endif DO 14 J=1,MFIT DO 13 K=1,MFIT COVAR(J,K)=0. 13 CONTINUE BETA(J)=0. 14 CONTINUE DO 18 I=1,NDATA CALL FUNCS(X(I),AFUNC,MA) YM=Y(I) IF(MFIT.LT.MA) THEN DO 15 J=MFIT+1,MA YM=YM-A(LISTA(J))*AFUNC(LISTA(J)) 15 CONTINUE ENDIF SIG2I=1./SIG(I)**2 DO 17 J=1,MFIT WT=AFUNC(LISTA(J))*SIG2I DO 16 K=1,J COVAR(J,K)=COVAR(J,K)+WT*AFUNC(LISTA(K)) 16 CONTINUE BETA(J)=BETA(J)+YM*WT 17 CONTINUE 18 CONTINUE IF (MFIT.GT.1) THEN DO 21 J=2,MFIT DO 19 K=1,J-1 COVAR(K,J)=COVAR(J,K) 19 CONTINUE 21 CONTINUE ENDIF CALL GAUSSJ(COVAR,MFIT,NCVM,BETA,1,1) DO 22 J=1,MFIT A(LISTA(J))=BETA(J) 22 CONTINUE CHISQ=0. DO 24 I=1,NDATA CALL FUNCS(X(I),AFUNC,MA) SUM=0. DO 23 J=1,MA SUM=SUM+A(J)*AFUNC(J) 23 CONTINUE CHISQ=CHISQ+((Y(I)-SUM)/SIG(I))**2 24 CONTINUE CALL COVSRT(COVAR,NCVM,MA,LISTA,MFIT) 25 RETURN END ************************************************************************ SUBROUTINE COVSRT(COVAR,NCVM,MA,LISTA,MFIT) *********************************************************************** CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C numerical recipes routine, Press et. al. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DOUBLE PRECISION COVAR(NCVM,NCVM),SWAP INTEGER*4 NCVM,MA,LISTA(MA),MFIT,I, * J DO 12 J=1,MA-1 DO 11 I=J+1,MA COVAR(I,J)=0. 11 CONTINUE 12 CONTINUE DO 14 I=1,MFIT-1 DO 13 J=I+1,MFIT IF(LISTA(J).GT.LISTA(I)) THEN COVAR(LISTA(J),LISTA(I))=COVAR(I,J) ELSE COVAR(LISTA(I),LISTA(J))=COVAR(I,J) ENDIF 13 CONTINUE 14 CONTINUE SWAP=COVAR(1,1) DO 15 J=1,MA COVAR(1,J)=COVAR(J,J) COVAR(J,J)=0. 15 CONTINUE COVAR(LISTA(1),LISTA(1))=SWAP DO 16 J=2,MFIT COVAR(LISTA(J),LISTA(J))=COVAR(1,J) 16 CONTINUE DO 18 J=2,MA DO 17 I=1,J-1 COVAR(I,J)=COVAR(J,I) 17 CONTINUE 18 CONTINUE RETURN END **************************************************************************** SUBROUTINE GAUSSJ(A,N,NP,B,M,MP) *************************************************************************** CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C numerical recipes routine, Press et. al. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC integer*4 NMAX PARAMETER (NMAX=50) INTEGER*4 IPIV(NMAX),INDXR(NMAX),INDXC(NMAX), * N,NP,M,MP,I,J,K,IROW,ICOL,L,LL DOUBLE PRECISION A(NP,NP),B(NP,NP),BIG,DUM, * PIVINV DO 11 J=1,N IPIV(J)=0 11 CONTINUE DO 22 I=1,N BIG=0. DO 13 J=1,N IF(IPIV(J).NE.1)THEN DO 12 K=1,N IF (IPIV(K).EQ.0) THEN IF (ABS(A(J,K)).GE.BIG)THEN BIG=ABS(A(J,K)) IROW=J ICOL=K ENDIF ELSE IF (IPIV(K).GT.1) THEN write(6,*) 'Singular matrix' goto 25 ENDIF 12 CONTINUE ENDIF 13 CONTINUE IPIV(ICOL)=IPIV(ICOL)+1 IF (IROW.NE.ICOL) THEN DO 14 L=1,N DUM=A(IROW,L) A(IROW,L)=A(ICOL,L) A(ICOL,L)=DUM 14 CONTINUE DO 15 L=1,M DUM=B(IROW,L) B(IROW,L)=B(ICOL,L) B(ICOL,L)=DUM 15 CONTINUE ENDIF INDXR(I)=IROW INDXC(I)=ICOL IF (A(ICOL,ICOL).EQ.0.) then write(6,*) 'Singular matrix.' goto 25 endif PIVINV=1./A(ICOL,ICOL) A(ICOL,ICOL)=1. DO 16 L=1,N A(ICOL,L)=A(ICOL,L)*PIVINV 16 CONTINUE DO 17 L=1,M B(ICOL,L)=B(ICOL,L)*PIVINV 17 CONTINUE DO 21 LL=1,N IF(LL.NE.ICOL)THEN DUM=A(LL,ICOL) A(LL,ICOL)=0. DO 18 L=1,N A(LL,L)=A(LL,L)-A(ICOL,L)*DUM 18 CONTINUE DO 19 L=1,M B(LL,L)=B(LL,L)-B(ICOL,L)*DUM 19 CONTINUE ENDIF 21 CONTINUE 22 CONTINUE DO 24 L=N,1,-1 IF(INDXR(L).NE.INDXC(L))THEN DO 23 K=1,N DUM=A(K,INDXR(L)) A(K,INDXR(L))=A(K,INDXC(L)) A(K,INDXC(L))=DUM 23 CONTINUE ENDIF 24 CONTINUE 25 RETURN END ******************************************************************************** Double Precision FUNCTION FRAC ( X1, X2, F) ******************************************************************************** C C Interpolates to find proper velocity width boundaries C Double Precision X1, X2, F C FRAC = (F-X1)/(X2-X1) c RETURN END