PROGRAM BetaCorner ! Program to calculate likelihood function for the two adjustable ! parameters (beta, and corner magnitude) of a tapered Gutenberg-Richter ! distribution [proposed by D. Vere-Jones; see David D. Jackson and ! Yan Y. Kagan, Testable earthquake forecasts for 1999, Seism. Res. Lett., ! volume 70, pages 393-403, 1999 for formal statement.]. ! Written in May of 1999 ! by Yan Y. Kagan, Institute of Geophysics and Planetary Physics, ! University of California, Los Angeles, CA 90095. ! Modified 2003.02 by Peter Bird, Department of Earth and Space ! Sciences, UCLA: ! -Fortran 66 fixed-format code converted to Fortran 90 free-format. ! -User-specified PARAMETERs collected together and annotated. ! (See the portion of code below which is outlined by # characters: ! ########################## ! FILE = '?' ! magnitude_FORMAT = '?' ! threshold_magnitude = ? ! tlab = '?' ! ########################## ! This is where you specify the name of the catalog input ! file, to be READ on unit 20. The earthquake catalog can ! be in any format as long as it includes a precise magnitude, ! and has one earthquake per line of the catalog file.) ! -Output files of text and numbers are given names: ! BetaCorner_03.txt ! BetaCorner_04.txt ! BetaCorner_06.txt (another copy of what the user sees) ! BetaCorner_11.txt ! -CALL Pause() added at end to be sure user has a chance ! to read the output from unit 6. (Some systems ! will close this window before it can be read.) ! -Graphics provided by my MODULE Adobe_Illustrator, not by GKS. ! (This output will appear as file BetaCorner.ai on unit 101. ! The supporting file AI7frame.ai must be available ! in the current directory, and will be read on unit 100.) ! -Unused SUBROUTINE Agutol was deleted. ! Additional notes below are by Yan Kagan: ! ----------------------------------------------------------------- ! 1999/11/24 realized that this calculation is exact, since ! expression for pdf does not have gamma function, whereas ! for gamma distribution the pdf formula has a normalizing ! coefficient with incomplete gamma function ! (Kagan, GJI, 106, 123-134, 1991, eqs. 6,7), DZGD2.FOR uses ! approximation for incomplete gamma function ! 2000/10/01 ! Added a moment formula for Mmax (Rick's) with bias corrected ! 2000/10/05 ! Corrected ML formula for truncated Pareto distribution (09/26) ! 2000/10/06 ! ML formula for truncated Pareto distribution can be used for ! unbiased estimate of Mc (formula by Pisarenko-Kijko): ! uncomment AMOMXC = AMOMX + AMOMX*(AMOMX**BETA - 1.0)/(BETA*SUMAL) ! Run first DZM.FOR, DZMM.FOR, or DZMM1.FOR program (CMT catalog) ! or TRANSFORM.FOR (Ocean transform earthquakes) ! Takes about 25 minutes for N=3500 super-threshold earthquakes on CYCLOP, ! but only 3-5 seconds for N=1920 super-threshold earthquakes on a 3 GHz Pentium-4 computer. ! --------------------------------------------------------------------------------------------- IMPLICIT REAL * 8 (a - h, o - z) PARAMETER (betau = 2.0d0 / 3.0d0) ! a priori estimate of beta, based on experience INTEGER, PARAMETER :: nxt = 91, & ! nxt = number of beta values to test & nyt = 91 ! nyt = number of corner magnitudes to test PARAMETER (nxttb = 31, nyttb = 16) ! nxttb is never used; nyttb is used only to limit the number of columns (j-values) of "table" printed INTEGER, PARAMETER :: catalog_limit = 30000 ! should exceed count of earthquakes in largest anticipated catalog DIMENSION amom(catalog_limit) ! list of moment ratios (with respect to threshold DIMENSION table (nxt, nyt) ! table of [first] absolute, [then] relative log-likelihoods (i = 1:nxt = beta index, j = 1:nyt = corner magnitude index) DIMENSION bline(nxt) ! list of values of beta corresponding to 1st subscript of "table" DIMENSION aline(nyt) ! list of values of LOG10(corner moment, in N m) corresponding to 2nd subscript of "table" DIMENSION alinm(nyt) ! list of corner magnitudes corresponding to 2nd subscript of "table" CHARACTER*50 magnitude_FORMAT ! text string describing how to get magnitude from each line of the catalog; e.g., "(53X, F5.2)" CHARACTER*50 tlab ! text label for plot produced by CALL Autoplot REAL :: figure_width_inches, figure_height_inches LOGICAL :: bitmap_wanted ! REAL*4 mag(8) !Engdahl catalog ! character*2 msc(8) !Engdahl catalog ! character*5 mdo(8), icat,asol,isol !Engdahl catalog !##################### USER: SELECT PARAMETER VALUES HERE !!! ###################################### ! Insert file name of seismic catalog in "FILE =" field. ! A seismic catalog file has one earthquake per line, and ! the only necessary information is the earthquake magnitude, ! which should always be in the same columns, and should be ! specified with 3 or more significant digits ! (e.g., " 7.13" is better than " 7.1", and " 7" is unacceptable. OPEN (UNIT = 20, STATUS = 'OLD', disp = 'KEEP', readonly, & & FILE = "merged_PB2002_SUB_all.cat") ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]ccb_pure.DAT') ! P. Bird catalogs ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]ctf_pure.DAT') ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]crb_pure.DAT') ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]osr_pure.DAT') ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]otf_pure.DAT') ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]ocb_pure.DAT') ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]sub_pure.DAT') ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]int_pure.DAT') ! & FILE= 'EQUAKEDKA0:[KAGAN]temp.DAT') ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]OTF_PURE_2-34MMPA.DAT') ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]OTF_PURE_35-66MMPA.DAT') ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]OTF_PURE_67-263MMPA.DAT') ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]MERG_CCB_ALL.DAT') ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]MERG_CTF_ALL.DAT') ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]MERG_OCB_ALL.DAT') ! & FILE= 'EQUAKEDKA0:[KAGAN.CAT]MERG_SUB_ALL.DAT') !------------------------------------------------------------------ ! "magnitude_FORMAT" is the Fortran FORMAT for reading the earthquake magnitude from the catalog file: magnitude_FORMAT = "(53X, F5.2)" ! P. Bird's .cat files !------------------------------------------------------------------ ! ! "threshold_magnitude": ! In recent papers, we call the threshold magnitude "m sub t". ! Earthquakes with magnitude <= threshold magnitude ! will be passed over, and not used in the calculation. Only ! earthquakes with magnitude > threshold_magnitude will be used. ! Generally, set the threshold magnitude to the smallest ! value for which the catalog is thought to be complete! ! It is also important to use a double-precision literal ! (e.g., "5.70D0") for the value of threshold_magnitude, ! because it will be compared to F5.2 values assigned to ! double-precision variable amag. If you use a single- ! precision literal (e.g., "5.70") then its value will be ! slightly less, and earthquakes with magnitude equal to ! threshold_magnitude will probably be admitted, causing a ! slightly different answer. Worse, results might even be ! compiler-dependent. ! Another way to prevent such threshold-precision problems ! is to specify threshold_magnitude as an in-between value ! that can never occur in the data set (e.g., "5.705D0"). ! threshold_magnitude = 5.40D0 ! CRB, OSR (from CMT) ! threshold_magnitude = 5.50D0 ! OTF (from CMT) ! threshold_magnitude = 5.70D0 ! CCB, CTF, OCB, SUB, INT (from CMT) threshold_magnitude = 7.10D0 ! merged (from Pacheco & Sykes [1992] + Ekstrom & Nettles [1997] + CMT (m>7)) !------------------------------------------------------------------ ! "tlab" is the text label (CHARACTER*50) for the output plot: ! tlab='CCB-pure 1977-2002.9.30, m>5.7' ! tlab='CTF-pure 1977-2002.9.30, m>5.7' ! tlab='CRB-pure 1977-2002.9.30, m>5.4' ! tlab='OSR-pure 1977-2002.9.30, m>5.4' ! tlab='OSR-pure-normal 1977-2002.9.30, m>5.4' ! tlab='OSR-pure-other 1977-2002.9.30, m>5.4' ! tlab='OTF-pure 3-39 mm/a 1977-2002.9.30, m>5.5' ! tlab='OTF-pure 40-68 mm/a 1977-2002.9.30, m>5.5' ! tlab='OTF-pure 69-263 mm/a 1977-2002.9.30, m>5.5' ! tlab='OCB-pure 1977-2002.9.30, m>5.7' ! tlab='SUB-pure 1977-2002.9.30, m>5.7' ! tlab='INT-pure 1977-2002.9.30, m>5.7' ! tlab='merged-CCB-all 1900-2002.9.30, m>7.1' ! tlab='merged-CTF-all 1900-2002.9.30, m>7.1' ! tlab='merged-CRB-all 1900-2002.9.30, m>7.1' ! tlab='merged-OTF 3-39 mm/a 1900-2002.9.30, m>7.1' ! tlab='merged-OCB-all 1900-2002.9.30, m>7.1' tlab='merged-SUB-all 1900-2002.9.30, m>7.1' !------------------------------------------------------------------ ! "betas" [for "beta starting value"] is the lowest value of beta to be tested. ! When in doubt, use a low value like 0.10D0. Zero and negative values are impossible. betas = 0.10D0 !simulation !------------------------------------------------------------------ ! "bmax" [for "beta maximum value"] is the largest value of beta to be tested. ! When in doubt, use 0.999D0, as beta >= 1 is theoretically impossible. bmax = 0.999D0 !simulation, Engdahl catalog !------------------------------------------------------------------ ! "ast" [for "a start"] is the initial value (when j==1 in table) of ! a == LOG10(threshold moment / trial corner moment) ! This will determine the lowest trial corner magnitude tested. ! That is, as ast increases above 1, ! this program looks further and further below the threshold magnitude ! in searching for the optimum corner magnitude; ! as ast decreases below 1, ! this program begins the search for corner magnitude ! progressively higher up from the threshold magnitude. ! Note that it take a change by 10.0**(+-1.5) to change the ! range of corner magnitudes considered by -+1.0 magnitude units. ast = 10.0**(-1.5) !Bird57_CCB_pure, Bird57_SUB_pure, Bird71_merged_SUB_all, Bird57_INT_pure ! ast = 10.0**(-1.0) !Bird57_CTF_pure, Bird54_CRB_pure, Bird55_OTF_pure, Bird57_OCB_pure ! ast = 10.0**(-0.5) ! ! ast = 1.0 !default; Bird71_merged_CCB_all, Bird71_merged_CTF_all, Bird71_merged_OTF_all_slow, Bird71_merged_OCB_all ! ast = 10.0**(0.5) !Bird71_merged_CRB_all, Bird54_OSR_pure_normal, Bird54_OSR_pure_other !------------------------------------------------------------------ ! "stepa" is the size of the step in LOG10(trial corner moment, in N m) ! when going from row j to row (j+1) in "table"; ! the corresponding step in corner magnitude is (2/3)*stepa ! stepa = 0.100D0 !OSR_pure_other; corner magnitude step 0.0666666... stepa = 0.050D0 !ALL except OSR_pure_other; corner magnitude step 0.033333... ! stepa = 0.0250D0 !corner magnitude step 0.016666... ! stepa = 0.0100D0 !corner magnitude step 0.006666... !######################## END OF PARAMETER SELECTION ############################################### OPEN (UNIT = 3, FILE = "BetaCorner_03.txt") ! unconditional; overwrites any existing file OPEN (UNIT = 4, FILE = "BetaCorner_04.txt") ! unconditional; overwrites any existing file OPEN (UNIT = 106, FILE = "BetaCorner_06.txt") ! unconditional; overwrites any existing file OPEN (UNIT = 11, FILE = "BetaCorner_11.txt") ! unconditional; overwrites any existing file !Variables derived from the user-selected parameters: stepb = (bmax - betas) / (nxt - 1) ! size of step in beta. NOTE that highest attainable precision in result is +- (stepb / 2). cutmap = -15.5 ! lower limit on values of "table" to be displayed. amcof = 1.5 thres = threshold_magnitude * 1.5 + 16.0 ! thres = LOG10(threshold moment, in dyne cm) amcut = threshold_magnitude amcut0 = amcof * 10 * threshold_magnitude ! amcut0 = (3/2) * 10 * threshold magnitude amomc = 10.0D0**(amcut0 / 10.0D0) ! amomc = 10**((3/2) threshold magnitude); ! this is "almost" the threshold moment, except it lacks ! the constant factor of 10**9 to 10**9.05 (N m, or SI units), ! or 10**16 to 10**16.05 (dyne cm, or cgs units). !Initialize counters before reading catalog: nn = 0 ! number of earthquakes in catalog nmom = 0 ! number of earthquakes exceeding threshold sumal = 0.0 ! REAL variable equal to nmom; number of earthquakes exceeding threshold amomx = 0.0 ! largest value of "amomt" in catalog; ! amomt = (current earthquake scalar moment)/(scalar moment of threshold magnitude) summom = 0.0 ! sum of: (earthquake moment)/(threshold moment) sumlog = 0.0 ! sum of: ln((earthquake moment)/(threshold moment)) summomd = 0.0 ! sum of: ((earthquake moment)/(threshold moment))**2 30 CONTINUE ! beginning of catalog-reading loop; will GO TO this line from below READ (20, magnitude_FORMAT, END = 40) amag nn = nn + 1 ! update count of earthquakes in catalog WRITE (11, 111) nn, amomt 111 FORMAT (I8, F11.3) IF (amag <= threshold_magnitude) GO TO 30 ! if earthquake is .LE. threshold, ignore it (and go read in another one) amomt = 10.0D0**((amag - amcut) * 1.5) ! amomt = (scalar moment of current catalog earthquake)/(scalar moment of threshold earthquake) ! ! NOTE: Various authors use slightly different values of the conversion constant ! C [of Kagan, 2003, PEPI, v. 135, p. 173-209] ! in converting from magnitude to moment, and vice versa. ! Here the earthquake moment is computed relative to the moment of the threshold ! event, which gets around the problem. All that is necessary is that the ! same magnitude scale should be used in defining threshold_magnitude ! as is used in the seismic catalog. ! The maximum-likelihood estimate of the corner magnitude output by this ! program will automatically be expressed using the same version of the magnitude scale. ! ! However, scalar moments computed within this program are according to ! one particular scale, which is based on an apparent definition in ! Hanks and Kanamori [1979]: ! m = (2/3) (log M - 9.05) ! which, inverted, becomes: ! M = 10**(9.05 + (3/2) * m)) ! (All formulas in this note assume that moment is in N m, according to the SI system.) ! Therefore, the maximum-likelihood corner moment output from this program ! **MAY NOT** be on the same scale as that assumed by the user. ! Fortunately, the maximum possible difference is small (12% in moment; 0.0333 in magnitude). nmom = nmom + 1 ! count one more earthquake exceeding threshold IF (nmom <= catalog_limit) THEN amom(nmom) = amomt ! add to list of saved moment ratios (to threshold moment) ELSE WRITE (6, "(' Increase catalog_limit (currently = ',I6,') and recompile.')") catalog_limit WRITE (106, "(' Increase catalog_limit (currently = ',I6,') and recompile.')") catalog_limit CALL Pause() STOP END IF IF (amomx.LT.amomt) amomx = amomt ! update saved value of largest amomt sumal = sumal + 1.0 ! increment REAL count of super-threshold earthquakes summom = summom + amomt ! increment sum of: (earthquake moment)/(threshold moment) summomd = summomd + amomt * amomt ! increment sum of: ((earthquake moment)/(threshold moment))**2 sumlog = sumlog + DLOG(amomt) ! increment sum of: ln((earthquake moment)/(threshold moment)) GO TO 30 ! loop to read another earthquake from catalog 40 CONTINUE ! exit point; will GO TO this line when any READ fails CLOSE (20) ! whole catalog has now been READ amean = summom / sumal ! mean ratio of [super-threshold] earthquake moment to threshold moment stdev = DSQRT(summomd / sumal - amean**2) ! standard deviation of [super-threshold] (earthquake moment)/(threshold moment) ammaxu = (summomd / sumal - 1.0) / (2.0 * (amean + betau - amean * betau)) bias = ( (2.0 + 3.0 * ammaxu * betau + (stdev**2 + amean**2) * & & (6.0 * ammaxu - 3.0 * ammaxu * betau - 2.0 * amean) ) * (1.0 - & & betau) ) / (4.0 * sumal * (betau + (1.0 - betau) * amean)**2) ammaxua = ammaxu + bias cmmaxua = 10.0**thres * ammaxua ammagu = (DLOG10(cmmaxua) - 16.0) / 1.5 WRITE (6, 233) amean, stdev, ammaxu, ammaxua, cmmaxua, ammagu WRITE (106, 233) amean, stdev, ammaxu, ammaxua, cmmaxua, ammagu 233 FORMAT (' amean, stdev, ammaxu, ammaxua, cmmaxua, ammagu =', / , 1P6G15.7) tabm = -1.0D20 ! highest log-likelihood in "table"; will be replaced amax = -1.0D20 betax = -1.0D20 tabmgr = -1.0D20 tabmgrt = -1.0D20 n = 0 ! count of entries in "table" that have been completed WRITE (6, 1) WRITE (106, 1) 1 FORMAT ('1') amomxc = amomx + amomx * (amomx**betau - 1.0) / (betau * sumal) WRITE (6, 168) amcut0, amomc, amomx, amomxc WRITE (106, 168) amcut0, amomc, amomx, amomxc 168 FORMAT (' AMCUT0, AMOMC, AMOMX, amomxc = ', 1P4G20.7) !Begin main computational loops in this program: DO 100 i = 1, nxt ! trial values of beta, increasing from betas to bmax: beta = betas + stepb * DFLOAT(i - 1) ! current trial value of beta bline(i) = beta ! store trial values of beta in list "bline" gr = sumal * DLOG(beta) - (1.0 + beta) * sumlog amomxc = amomx ! amomxc = amomx + amomx*(amomx**beta - 1.0)/(beta*sumal) !truncate grt = sumal * (beta * DLOG(amomxc) - DLOG(amomxc**beta - 1.0) + & & DLOG(beta) ) - (1.0 + beta) * sumlog IF (gr > tabmgr) tabmgr = gr IF (grt > tabmgrt) tabmgrt = grt DO 110 j = 1, nyt ! trial values of corner magnitude [or log corner moment]. a = ast * 10.0**(DFLOAT(1 - j) * stepa) ! when j == 1, initial value is ast ["a start"] aline(j) = -DLOG10(a) + 9.05 + DLOG10(amomc) ! LOG10(trial corner moment, in N m), increasing with j alinm(j) = (aline(j) - 9.05) / amcof ! trial corner magnitude, increasing with j ! GAM = DGAMMA(1.0D0-BETA) ! F1 = SUMAL*DLOG(BETA) !C F2 = SUMAL*GAM*A**BETA ! F2 = SUMAL*DLOG(1.0D0 - A**BETA*DEXP(A)*(GAM - ! 1 (A**(1.0D0-BETA)/(1.0D0-BETA)) + (A**(2.0D0-BETA)/(2.0D0-BETA)) ! 2 (A**(3.0D0-BETA)/(2.0D0*(3.0D0-BETA))))) ! F3 = A*SUMMOM ! F4 = (1.0D0 + BETA)*SUMLOG ! F5 = SUMAL*A !c TABL = F1 + F2 - F3 - F4 ! TABL = F1 - F2 - F3 - F4 + F5 n = n + 1 f2a = 0.0 DO 10 imom = 1, nmom f2a = f2a + DLOG(beta / amom(imom) + a) 10 CONTINUE f3 = a * summom f4b = beta * sumlog f5 = sumal * a tabl = f2a - f3 - f4b + f5 IF (MOD(n, 250) == 1) THEN WRITE (6, 200) i, j, a, beta, f2a, f3, f4b, f5, tabl WRITE (106, 200) i, j, a, beta, f2a, f3, f4b, f5, tabl 200 FORMAT (' I,J,A,BETA,F2A,F3,F4B,F5,TABL= ', /, 2I3, 9G11.4) END IF table (i, j) = tabl ! store absolute log-likelihood in "table"(i = 1:nxt for beta, j = 1:nyt for corner moment) IF (tabl > tabm) THEN tabm = tabl ! tabm = latest estimate of highest log-likelihood in "table" amax = a ! amax = latest estimate of "a" for highest log-likelihood point in "table" betax = beta ! bmax = latest estimate of beta for highest log-likelihood point in "table" imax = i ! imax = latest estimate of 1st subscript at highest log-likelihood point in "table" jmax = j ! jmax = latest estimate of 2nd subscript at highest log-likelihood point in "table" END IF 110 CONTINUE 100 CONTINUE WRITE (6, 210) imax, jmax, nmom, tabm, tabmgr, tabm - tabmgr, & & tabmgrt, tabm - tabmgrt WRITE (106, 210) imax, jmax, nmom, tabm, tabmgr, tabm - tabmgr, & & tabmgrt, tabm - tabmgrt 210 FORMAT ('0 IMAX, JMAX, nmom, TABM, TABMGR, Diffgr, TABMGRT, ', & & 'Diffgrt = ', /, 3I7, 5G19.9, //) temp2 = amomc / amax amagmax = DLOG10(temp2) / 1.5 WRITE (3, 112) nn, nmom, amcut, betax, amagmax 112 FORMAT (2i8, 3f12.5) WRITE (6, 215) amax, 1.0 / amax, betax, temp2, amagmax WRITE (106, 215) amax, 1.0 / amax, betax, temp2, amagmax 215 FORMAT ('0 AMAX, 1.0/AMAX, betamax, mommax, magmax = ', & & /, 1p7G15.7) WRITE (6, 1) WRITE (106, 1) WRITE (6, 240) bline WRITE (106, 240) bline WRITE (6, 220) WRITE (106, 220) WRITE (6, 240) aline WRITE (106, 240) aline WRITE (6, 220) WRITE (106, 220) WRITE (6, 240) alinm WRITE (106, 240) alinm 240 FORMAT (' ', 12X, 10F12.4) WRITE (6, 220) WRITE (106, 220) WRITE (6, 220) WRITE (106, 220) WRITE (6, 220) WRITE (106, 220) !Calibrate array "table" so that its maximum log-likelihood becomes +3.0D0: DO 120 i = 1, nxt DO 130 j = 1, nyt table (i, j) = table (i, j) - tabm + 3.0D0 IF (table (i, j).LT. - 999.99) table (i, j) = -999.99 130 CONTINUE 120 CONTINUE !Note that "table" now contains relative, not absolute, log-likelihoods. !Output preview of "table": ! successive rows have increasing i, or beta; ! successive columns have increasing j, or corner magnitude. ! Note that not all values of j/corner magnitude are output ! to the screen, but all go to unit 106. WRITE (6, 310) (aline(i), i = 1, nyttb) WRITE (106, 311) (aline(i), i = 1, nyt) WRITE (6, 310) (alinm(i), i = 1, nyttb) WRITE (106, 311) (alinm(i), i = 1, nyt) WRITE (106, *) 310 FORMAT (' ', 9X, 16F7.2, /) 311 FORMAT (' ', 9X, 99F7.2) ! 99 may be more than needed; no harm done DO 160 i = 1, nxt beta = betas + stepb * DFLOAT(i - 1) WRITE (6, 220) beta, (table(i, j), j = 1, nyttb) 220 FORMAT (' ', F7.3, 2X, 16F7.2) WRITE (106, 221) beta, (table(i, j), j = 1, nyt) 221 FORMAT (' ', F7.3, 2X, 99F7.2) ! 99 may be more than needed; no harm done 160 CONTINUE WRITE (6, 1) WRITE (106, 1) xmin = betas xmax = betas + stepb * DFLOAT(nxt - 1) ymin = aline(1) ymax = aline(nyt) ! figure_width_inches = 4.00 ! small, to fit Yan's notebooks figure_width_inches = 2.30 ! very small, to fit 6 on one 8.5*11" page figure_height_inches = figure_width_inches ! square format bitmap_wanted = .TRUE. ! IF(bitmap_wanted), figure includes color representation of log-likelihood. CALL Autoplot (xmin, xmax, ymin, ymax, tabm, tabmgr, tabmgrt, & ! original Kagan parameters & betax, amagmax, ammagu, amcut, nmom, & ! original Kagan parameters & table, nxt, nyt, & ! extra parameters added by Bird & tlab, & ! extra parameters added by Bird & figure_width_inches, figure_height_inches, & ! extra parameters added by Bird & bitmap_wanted) ! extra parameters added by Bird !PB CALL Conrec (table, nxt, nxt, nyt, 0., 0., 3., -1, 0, -992) !shall !PB CALL Frame WRITE (6, 133) nmom, sumal, sumlog, summom WRITE (106, 133) nmom, sumal, sumlog, summom 133 FORMAT (' N, SUMAL, SUMLOG, SUMMOM =', I5, F12.3, 2F13.3) CALL Pause() STOP CONTAINS SUBROUTINE Autoplot (xmin, xmax, ymin, ymax, tabm, tabmgr, & ! original Kagan parameters & tabmgrt, betax, amagmax, ammagu, amcut, nmom, & ! original Kagan parameters & table, nxt, nyt, & ! extra parameters added by Bird & tlab, & ! extra parameters added by Bird & figure_width_inches, figure_height_inches, & ! extra parameters added by Bird & bitmap_wanted) ! extra parameters added by Bird ! Use MODULE Adobe_Illustrator support routines !(from Peter Bird's file Adobe_Illustrator.f90) ! to create a one-page Adobe Illustrator graphic file (BetaCorner.ai) ! of contours of the log-likelihood values in "table", ! annotated with information about the axes and maximum-likelihood result. USE Adobe_Illustrator DIMENSION xdra(4), ydra(4) CHARACTER*50 tlab CHARACTER likmax*18 CHARACTER numb*18 CHARACTER gr*18 CHARACTER grt*18 CHARACTER magx*18 CHARACTER betx*18 CHARACTER magu*18 REAL*8 tabm, tabmgr, tabmgrt, betax, amagmax, ammagu, amcut INTEGER, INTENT(IN) :: nxt, nyt REAL*8, DIMENSION(nxt, nyt), INTENT(IN) :: table ! table of log-likelihoods (i = 1:nxt = beta index, j = 1:nyt = corner magnitude index) REAL, INTENT(IN) :: figure_width_inches, figure_height_inches LOGICAL, INTENT(IN) :: bitmap_wanted ! should figure be in color? LOGICAL :: in_ok, out_ok REAL, PARAMETER :: text_border_inches = 0.4 ! between outside edge of contour-map window and inside edge of blank margins REAL, PARAMETER :: tick_length_inches = 0.06 REAL :: top_margin_points, left_margin_points, right_margin_points, bottom_margin_points REAL :: x1_points, x2_points, y1_points, y2_points REAL :: bitmap_color_lowvalue, bitmap_color_highvalue, brightness, t CHARACTER(LEN=3) :: c3 CHARACTER(LEN=5) :: c5 CHARACTER(LEN=3), DIMENSION(:,:), ALLOCATABLE :: bitmap INTEGER, PARAMETER :: maximum_bits = 1000 ! should be plenty if grids have < 300 rows/columns INTEGER, PARAMETER :: maximum_labels = 100 ! should be plenty, as only ~6 are expected INTEGER :: bits_left, cell_intersections, & & far_end, grab_end, & & i, i_next, i_start, j, & & n_bits,n_contour, n_labels, n_wide, n_x_cells, n_y_cells, & & pixel_rows, pixel_columns, & & x_index_high, x_index_low, y_index_high, y_index_low LOGICAL :: pen_down LOGICAL, DIMENSION(4) :: cell_free_end LOGICAL, DIMENSION(maximum_bits) :: used_bit LOGICAL, DIMENSION(2, maximum_bits) :: free_end_12 REAL :: beta, beta_tick_interval, & & corner_magnitude_tick_interval, criterion, criterion_bestYet, current_x_points, current_y_points, & & dx_points, dx_prime, dy_points, dy_prime, & & fraction, & & gradient, & & label_angle_radians, label_x_points, label_y_points, left_value, log10_M_c_tick_interval, loglike, & & maximum_corner_magnitude, minimum_corner_magnitude, minimum_distance_points, & & right_value, & & tick_length_points, t11, t12, t21, t22, tx, txp, tx1, ty, typ, ty1, ty2, ty3, ty4, ty5, ty6, & & x_high_points, x_index_fraction, x_low_points, x1_box_points, x2_box_points, xp, xp1, xp2, & & y_high_points, y_index_fraction, y_low_points, y1_box_points, y2_box_points, yp, yp1, yp2 REAL, DIMENSION(2, 4) :: cell_xy_points, xy_4_corner_points REAL, DIMENSION(2, maximum_labels) :: label_xy_points REAL, DIMENSION(maximum_bits) :: gradient_x, gradient_y REAL, DIMENSION(2, 2, maximum_bits) :: bit_xy_12_points !Initialize Adobe_Illustrator by choosing paper size, margin size, and contour-map window placement: IF (figure_width_inches > figure_height_inches) THEN ! use landscape format CALL Select_Paper (paper_width_points = 11.00 * 72, paper_height_points = 8.25 * 72) ELSE ! use portrait format CALL Select_Paper (paper_width_points = 8.25 * 72, paper_height_points = 11.00 * 72) END IF CALL Set_Background(black = .FALSE.) left_margin_points = MAX(0.0, (ai_paper_width_points - (figure_width_inches * 72) - 2 * (text_border_inches * 72)) / 2.0) right_margin_points = left_margin_points top_margin_points = MAX(0.0, (ai_paper_height_points - (figure_height_inches * 72) - 2 * (text_border_inches * 72)) / 2.0) bottom_margin_points = top_margin_points CALL Define_Margins (top_margin_points, & & left_margin_points, right_margin_points, & & bottom_margin_points) CALL Begin_Page (model_ai_filename = "AI7frame.ai", in_ok = in_ok, & & new_ai_filename = "BetaCorner.ai", out_ok = out_ok, & & using_color = bitmap_wanted, & & plan_toptitles = .FALSE., & & plan_rightlegend = .FALSE., & & plan_bottomlegend = .FALSE.) IF (.NOT.in_ok) THEN WRITE (*, "(' ERROR: Required input file AI7frame.ai could not be found (in current folder).')") CALL Pause() STOP END IF IF (.NOT.out_ok) THEN WRITE (*, "(' ERROR: Output file BetaCorner.ai already exists. Please delete, move, or rename it.')") CALL Pause() STOP END IF x1_points = left_margin_points + (text_border_inches * 72) x2_points = ai_paper_width_points - right_margin_points - (text_border_inches * 72) y1_points = bottom_margin_points + (text_border_inches * 72) y2_points = ai_paper_height_points - top_margin_points - (text_border_inches * 72) CALL Set_Window (x1_points, x2_points, y1_points, y2_points) IF (bitmap_wanted) THEN pixel_rows = NINT(3 * (x2_points - x1_points) + 1.) ! should be plenty pixel_columns = pixel_rows ! since contour-map window is square ALLOCATE( bitmap(pixel_rows, pixel_columns) ) bitmap_color_highvalue = +3.0 bitmap_color_lowvalue = 0.0 brightness = 1.3 ! use 1.0 for full saturation, >1.0 to mix in white, <1.0 to mix in black DO i = 1, pixel_rows ! note that these pixel rows count down from the top! y_index_high = Int_Above(1 + ((pixel_rows - i) * (nyt - 1.)) / (pixel_rows - 1.)) y_index_high = MIN(y_index_high, nyt) y_index_low = MAX((y_index_high - 1), 1) y_index_fraction = 1 + ((pixel_rows - i) * (nyt - 1.)) / (pixel_rows - 1.) - y_index_low ! increasing (0.0-1.0) between y_index_low and y_index_high DO j = 1, pixel_columns x_index_high = Int_Above(1 + (j - 1.) * (nxt - 1.) / (pixel_columns - 1.)) x_index_high = MIN(x_index_high, nxt) x_index_low = MAX((x_index_high - 1), 1) x_index_fraction = 1 + (j - 1.) * (nxt - 1.) / (pixel_columns - 1.) - x_index_low left_value = y_index_fraction * table(x_index_low , y_index_high) + (1. - y_index_fraction) * table(x_index_low , y_index_low) right_value = y_index_fraction * table(x_index_high, y_index_high) + (1. - y_index_fraction) * table(x_index_high, y_index_low) value = x_index_fraction * right_value + (1. - x_index_fraction) * left_value IF (value >= bitmap_color_lowvalue) THEN t = (value - bitmap_color_lowvalue) / (bitmap_color_highvalue - bitmap_color_lowvalue) c3 = RGB_Munsell(warmth = t, brightness = brightness) bitmap(i, j) = c3 ELSE ! 10%-gray bitmap(i, j) = CHAR(230)//CHAR(230)//CHAR(230) END IF END DO END DO CALL Bitmap_on_L1 (bitmap, x1_points, x2_points, y1_points, y2_points) DEALLOCATE ( bitmap ) END IF ! (bitmap_wanted) !Draw contours of log-likelihood in the contour-map window. ! Note: The topology of contours is quite complex in the general case: ! contours may expand into black areas if several grid values are exactly ! the contour level, and an abitrary (even) number of contours may emanate ! from this patch, etc., etc. ! This code ignores such complexities, and assumes that it is ! quite rare for grid points to have values exactly on a contour. ! IF they do, minor topological errors may result. CALL Set_Stroke_Color("foreground") n_labels = 0 label_xy_points = 0.0 ! whole array, to make debugging easier DO n_contour = 0, -15, -3 loglike = REAL(n_contour) IF (loglike == 0.0) THEN ! use solid line CALL Set_Line_Style (width_points = 2.0, dashed = .FALSE.) ELSE ! use dashed line CALL Set_Line_Style (width_points = 1.0, dashed = .TRUE., on_points = 6.0, off_points = 3.0) END IF !Build a list of countour bits: n_bits = 0 bit_xy_12_points = 0.0 ! whole array; makes debugging easier free_end_12 = .FALSE. ! whole array; makes debugging easier gradient_x = 0.0 ! whole array; makes debugging easier gradient_y = 0.0 ! whole array; makes debugging easier n_x_cells = nxt - 1 n_y_cells = nyt - 1 ! a "cell" is a rectangle outlined by 4 adjacent grid points DO i = 1, n_x_cells x_index_low = i x_index_high = i+1 x_low_points = x1_points + (x2_points - x1_points) * (x_index_low - 1.) / (nxt - 1.) x_high_points = x1_points + (x2_points - x1_points) * (x_index_high - 1.) / (nxt - 1.) DO j = 1, n_y_cells y_index_low = j y_index_high = j+1 y_low_points = y1_points + (y2_points - y1_points) * (y_index_low - 1.) / (nyt - 1.) y_high_points = y1_points + (y2_points - y1_points) * (y_index_high - 1.) / (nyt - 1.) t11 = table(x_index_low, y_index_low) t12 = table(x_index_low, y_index_high) t21 = table(x_index_high, y_index_low) t22 = table(x_index_high, y_index_high) IF ((MAX(t11, t12, t21, t22) > loglike).AND.(MIN(t11, t12, t21, t22) < loglike)) THEN !One or more contour bits should be found in this cell; start building list of control points: cell_intersections = 0 cell_xy_points = 0.0 ! whole (tiny) array cell_free_end = .FALSE. ! whole (tiny) array !bottom (low-y) side of cell: IF (((t11 - loglike) * (t21 - loglike)) < 0.0) THEN cell_intersections = cell_intersections + 1 fraction = (loglike - t11) / (t21 - t11) cell_xy_points(1, cell_intersections) = x_low_points + (x_high_points - x_low_points) * fraction cell_xy_points(2, cell_intersections) = y_low_points cell_free_end(cell_intersections) = (y_index_low == 1) END IF !right (high-x) side of cell: IF (((t21 - loglike) * (t22 - loglike)) < 0.0) THEN cell_intersections = cell_intersections + 1 fraction = (loglike - t21) / (t22 - t21) cell_xy_points(1, cell_intersections) = x_high_points cell_xy_points(2, cell_intersections) = y_low_points + (y_high_points - y_low_points) * fraction cell_free_end(cell_intersections) = (x_index_high == nxt) END IF !top (high-y) side of cell: IF (((t22 - loglike) * (t12 - loglike)) < 0.0) THEN cell_intersections = cell_intersections + 1 fraction = (loglike - t22) / (t12 - t22) cell_xy_points(1, cell_intersections) = x_high_points + (x_low_points - x_high_points) * fraction cell_xy_points(2, cell_intersections) = y_high_points cell_free_end(cell_intersections) = (y_index_high == nyt) END IF !left (low-x) side of cell: IF (((t12 - loglike) * (t11 - loglike)) < 0.0) THEN cell_intersections = cell_intersections + 1 fraction = (loglike - t12) / (t11 - t12) cell_xy_points(1, cell_intersections) = x_low_points cell_xy_points(2, cell_intersections) = y_high_points + (y_low_points - y_high_points) * fraction cell_free_end(cell_intersections) = (x_index_low == 1) END IF IF (MOD(cell_intersections, 2) == 1) cell_intersections = cell_intersections - 1 IF (cell_intersections == 2) THEN !typical case: connect the dots to make a contour bit IF (n_bits < maximum_bits) THEN n_bits = n_bits + 1 bit_xy_12_points(1, 1, n_bits) = cell_xy_points(1, 1) bit_xy_12_points(2, 1, n_bits) = cell_xy_points(2, 1) bit_xy_12_points(1, 2, n_bits) = cell_xy_points(1, 2) bit_xy_12_points(2, 2, n_bits) = cell_xy_points(2, 2) free_end_12(1, n_bits) = cell_free_end(1) free_end_12(2, n_bits) = cell_free_end(2) gradient_x(n_bits) = t21 - t11 gradient_y(n_bits) = t12 - t11 ELSE ! table has filled up WRITE (*, "(' ERROR: Increase INTEGER, PARAMETER :: maximum_bits = ',I6)") maximum_bits CALL Pause() STOP END IF END IF IF (cell_intersections == 4) THEN !This cell has 4 boundary points. !There is no way to know the "proper" topology of this saddle points, !so we arbitrarily connect the boundary points #1-#2, #3-#4: IF (n_bits < maximum_bits) THEN n_bits = n_bits + 1 bit_xy_12_points(1, 1, n_bits) = cell_xy_points(1, 3) bit_xy_12_points(2, 1, n_bits) = cell_xy_points(2, 3) bit_xy_12_points(1, 2, n_bits) = cell_xy_points(1, 4) bit_xy_12_points(2, 2, n_bits) = cell_xy_points(2, 4) free_end_12(1, n_bits) = cell_free_end(3) free_end_12(2, n_bits) = cell_free_end(4) gradient_x(n_bits) = t21 - t11 gradient_y(n_bits) = t12 - t11 ELSE ! table has filled up WRITE (*, "(' ERROR: Increase INTEGER, PARAMETER :: maximum_bits = ',I6)") maximum_bits CALL Pause() STOP END IF ! room in bits table, or not END IF ! 4 boundary points for this cell END IF ! this cell has any boundary points END DO ! j = 1, n_y_cells END DO ! i = 1, n_x_cells !Indefinate loop to link together contour bits into multi-bit segments used_bit = .FALSE. ! whole array bits_left = n_bits !Each pass through the following loop should set down the pen, draw a line, and then lift the pen: linking: DO pen_down = .FALSE. IF (bits_left <= 0) EXIT linking !Find a bit to start the line with: i_start = 0 ! don't have a starting bit yet !Try to find a starting bit with a free end: find_free: DO i = 1, n_bits IF (.NOT.used_bit(i)) THEN IF (free_end_12(1, i).OR.free_end_12(2, i)) THEN i_start = i EXIT find_free END IF END IF END DO find_free IF (i_start == 0) THEN !no unused bit with a free end was found; take first unused bit find_any: DO i = 1, n_bits IF (.NOT.used_bit(i)) THEN i_start = i EXIT find_any END IF END DO find_any END IF IF (free_end_12(1, i_start)) THEN ! begin line at end 1 CALL New_L12_Path (level = 2, x_points = bit_xy_12_points(1, 1, i_start), y_points = bit_xy_12_points(2, 1, i_start)) pen_down = .TRUE. CALL Line_To_L12 (x_points = bit_xy_12_points(1, 2, i_start), y_points = bit_xy_12_points(2, 2, i_start)) used_bit(i_start) = .TRUE. bits_left = bits_left - 1 current_x_points = bit_xy_12_points(1, 2, i_start) current_y_points = bit_xy_12_points(2, 2, i_start) ELSE ! start with end 2, whether free or not! CALL New_L12_Path (level = 2, x_points = bit_xy_12_points(1, 2, i_start), y_points = bit_xy_12_points(2, 2, i_start)) pen_down = .TRUE. CALL Line_To_L12 (x_points = bit_xy_12_points(1, 1, i_start), y_points = bit_xy_12_points(2, 1, i_start)) used_bit(i_start) = .TRUE. bits_left = bits_left - 1 current_x_points = bit_xy_12_points(1, 1, i_start) current_y_points = bit_xy_12_points(2, 1, i_start) END IF !initialize search for best label position (and orientation) on this contour: gradient = SQRT(gradient_x(i_start)**2 + gradient_y(i_start)**2) label_x_points = (bit_xy_12_points(1, 1, i_start) + bit_xy_12_points(1, 2, i_start)) / 2 label_y_points = (bit_xy_12_points(2, 1, i_start) + bit_xy_12_points(2, 2, i_start)) / 2 minimum_distance_points = MIN((label_x_points - x1_points), & & (x2_points - label_x_points), & & (label_y_points - y1_points), & & (y2_points - label_y_points), 36.0) DO j = 1, n_labels minimum_distance_points = MIN(minimum_distance_points, & & SQRT((label_x_points - label_xy_points(1, j))**2 + & & (label_y_points - label_xy_points(2, j))**2)) END DO criterion = minimum_distance_points / gradient criterion_bestYet = criterion label_angle_radians = ATAN(gradient_y(i_start) / gradient_x(i_start)) - 1.5708 IF (COS(label_angle_radians) < 0.0) label_angle_radians = label_angle_radians + 3.14159 add_bits: DO ! indefinite loop to add bits to existing line, if possible IF (bits_left <= 0) EXIT add_bits i_next = 0 ! don't have the next bit, yet find_connector: DO i = 1, n_bits IF (.NOT.used_bit(i)) THEN IF ((bit_xy_12_points(1, 1, i) == current_x_points) .AND.(bit_xy_12_points(2, 1, i) == current_y_points)) THEN i_next = i grab_end = 1 far_end = 2 EXIT find_connector END IF IF ((bit_xy_12_points(1, 2, i) == current_x_points) .AND.(bit_xy_12_points(2, 2, i) == current_y_points)) THEN i_next = i grab_end = 2 far_end = 1 EXIT find_connector END IF END IF END DO find_connector IF (i_next == 0) EXIT add_bits CALL Line_To_L12 (x_points = bit_xy_12_points(1, far_end, i_next), y_points = bit_xy_12_points(2, far_end, i_next)) used_bit(i_next) = .TRUE. bits_left = bits_left - 1 current_x_points = bit_xy_12_points(1, far_end, i_next) current_y_points = bit_xy_12_points(2, far_end, i_next) !consider whether this contour bit provides a better location for the contour label? gradient = SQRT(gradient_x(i_next)**2 + gradient_y(i_next)**2) txp = (bit_xy_12_points(1, 1, i_next) + bit_xy_12_points(1, 2, i_next)) / 2 typ = (bit_xy_12_points(2, 1, i_next) + bit_xy_12_points(2, 2, i_next)) / 2 minimum_distance_points = MIN((txp - x1_points), & & (x2_points - txp), & & (typ - y1_points), & & (y2_points - typ), 36.0) DO j = 1, n_labels minimum_distance_points = MIN(minimum_distance_points, & & SQRT((txp - label_xy_points(1, j))**2 + & & (typ - label_xy_points(2, j))**2)) END DO criterion = minimum_distance_points / gradient IF (criterion > criterion_bestYet) THEN criterion_bestYet = criterion label_x_points = txp label_y_points = typ label_angle_radians = ATAN(gradient_y(i_next) / gradient_x(i_next)) - 1.5708 IF (COS(label_angle_radians) < 0.0) label_angle_radians = label_angle_radians + 3.14159 END IF END DO add_bits IF (pen_down) CALL End_L12_Path (close = .FALSE., stroke = .TRUE., fill = .FALSE.) !label this contour: n_labels = n_labels + 1 IF (n_labels <= maximum_labels) THEN label_xy_points(1, n_labels) = label_x_points label_xy_points(2, n_labels) = label_y_points ELSE WRITE (*, "(' ERROR: Increase INTEGER, PARAMETER :: maximum_labels = ',I6)") maximum_labels CALL Pause() STOP END IF CALL Begin_Group() WRITE (c3, "(I3)") n_contour n_wide = LEN_TRIM(ADJUSTL(c3)) !create white box, first in upright position xp1 = label_x_points - (n_wide / 2.0) * 5 - 2 xp2 = label_x_points + (n_wide / 2.0) * 5 + 2 yp1 = label_y_points - 4 yp2 = label_y_points + 4 xy_4_corner_points(1, 1) = xp1; xy_4_corner_points(2, 1) = yp1 xy_4_corner_points(1, 2) = xp2; xy_4_corner_points(2, 2) = yp1 xy_4_corner_points(1, 3) = xp2; xy_4_corner_points(2, 3) = yp2 xy_4_corner_points(1, 4) = xp1; xy_4_corner_points(2, 4) = yp2 !now rotate white box: DO j = 1, 4 dx_points = xy_4_corner_points(1, j) - label_x_points dy_points = xy_4_corner_points(2, j) - label_y_points dx_prime = dx_points * COS(label_angle_radians) - dy_points * SIN(label_angle_radians) dy_prime = dx_points * SIN(label_angle_radians) + dy_points * COS(label_angle_radians) xy_4_corner_points(1, j) = label_x_points + dx_prime xy_4_corner_points(2, j) = label_y_points + dy_prime END DO !finally, plot white box: CALL Set_Fill_or_Pattern (use_pattern = .FALSE., color_name = "background") CALL New_L12_Path (level = 1, x_points = xy_4_corner_points(1, 1), y_points = xy_4_corner_points(2, 1)) CALL Line_To_L12 (x_points = xy_4_corner_points(1, 2), y_points = xy_4_corner_points(2,2)) CALL Line_To_L12 (x_points = xy_4_corner_points(1, 3), y_points = xy_4_corner_points(2,3)) CALL Line_To_L12 (x_points = xy_4_corner_points(1, 4), y_points = xy_4_corner_points(2,4)) CALL Line_To_L12 (x_points = xy_4_corner_points(1, 1), y_points = xy_4_corner_points(2,1)) CALL End_L12_Path (close = .TRUE., stroke = .FALSE., fill = .TRUE.) !put number in front of white box, for legibility: CALL Set_Fill_or_Pattern (use_pattern = .FALSE., color_name = "foreground") CALL L12_Text (level = 2, x_points = label_x_points, y_points = label_y_points, & & angle_radians = label_angle_radians, & & font_points = 10, lr_fraction = 0.5, ud_fraction = +0.4, & & text = TRIM(ADJUSTL(c3))) CALL End_Group() END DO linking END DO ! n_countour = 0, -15, -3 !Mark highest value(s) in grid: tick_length_points = tick_length_inches * 72 CALL Set_Line_Style (width_points = 1.0, dashed = .FALSE.) CALL Set_Stroke_Color("foreground") DO i = 1, nxt DO j = 1, nyt IF (table(i, j) == 3.00) THEN xp = x1_points + (x2_points - x1_points) * (i - 1.) / (nxt - 1.) yp = y1_points + (y2_points - y1_points) * (j - 1.) / (nyt - 1.) CALL Begin_Group() CALL New_L12_Path (level = 1, x_points = xp - tick_length_points, y_points = yp) CALL Line_To_L12 (x_points = xp + tick_length_points, y_points = yp) CALL End_L12_Path (close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL New_L12_Path (level = 1, x_points = xp, y_points = yp + tick_length_points) CALL Line_To_L12 (x_points = xp, y_points = yp - tick_length_points) CALL End_L12_Path (close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL End_Group() CALL L12_Text (level = 1, x_points = xp + tick_length_points, y_points = yp, angle_radians = 0.0, & & font_points = 10, lr_fraction = -0.2, ud_fraction = +0.4, & & text = "3") END IF END DO END DO !Put black border around contour-map window: CALL Set_Line_Style (width_points = 2.0, dashed = .FALSE.) CALL Set_Stroke_Color("foreground") CALL New_L12_Path (level = 1, x_points = x1_points, y_points = y1_points) CALL Line_To_L12 (x_points = x2_points, y_points = y1_points) CALL Line_To_L12 (x_points = x2_points, y_points = y2_points) CALL Line_To_L12 (x_points = x1_points, y_points = y2_points) CALL Line_To_L12 (x_points = x1_points, y_points = y1_points) CALL End_L12_Path (close = .TRUE., stroke = .TRUE., fill = .FALSE.) !Put tick marks and labels along x = beta axis: beta_tick_interval = 0.1 x_index_low = Int_Above(REAL(xmin) / beta_tick_interval) x_index_high = Int_Below(REAL(xmax) / beta_tick_interval) CALL Begin_Group() ! of tick marks CALL Set_Line_Style (width_points = 1.0, dashed = .FALSE.) CALL Set_Stroke_Color("foreground") DO i = x_index_low, x_index_high beta = i * beta_tick_interval xp = x1_points + (x2_points - x1_points) * (beta - xmin) / (xmax - xmin) CALL New_L12_Path (level = 1, x_points = xp, y_points = y1_points) CALL Line_To_L12 (x_points = xp, y_points = y1_points + tick_length_points) CALL End_L12_Path (close = .FALSE., stroke = .TRUE., fill = .FALSE.) CALL New_L12_Path (level = 1, x_points = xp, y_points = y2_points) CALL Line_To_L12 (x_points = xp, y_points = y2_points - tick_length_points) CALL End_L12_Path (close = .FALSE., stroke = .TRUE., fill = .FALSE.) END DO CALL End_Group() CALL Begin_Group() ! of text labels CALL Set_Fill_or_Pattern (use_pattern = .FALSE., color_name = "foreground") DO i = x_index_low, x_index_high beta = i * beta_tick_interval WRITE (c5, "(F5.2)") beta xp = x1_points + (x2_points - x1_points) * (beta - xmin) / (xmax - xmin) CALL L12_Text (level = 1, x_points = xp, y_points = y1_points, angle_radians = 0.0, & & font_points = 10, lr_fraction = 0.5, ud_fraction = +1.0, & & text = ADJUSTL(c5)) END DO CALL End_Group() CALL L12_Text (level = 1, x_points = ai_window_xc_points, y_points = y1_points, angle_radians = 0.0, & & font_points = 11, lr_fraction = 0.5, ud_fraction = +2.0, & & text = "Beta") !Put tick marks and labels along the left vertical axis: Log10 (Mc / 1 N m) log10_M_c_tick_interval = 1.0 y_index_low = Int_Above(REAL(ymin) / log10_M_c_tick_interval) y_index_high = Int_Below(REAL(ymax) / log10_M_c_tick_interval) CALL Begin_Group() ! of tick marks CALL Set_Line_Style (width_points = 1.0, dashed = .FALSE.) CALL Set_Stroke_Color("foreground") DO j = y_index_low, y_index_high t = j * log10_M_c_tick_interval yp = y1_points + (y2_points - y1_points) * (t - ymin) / (ymax - ymin) CALL New_L12_Path (level = 1, x_points = x1_points, y_points = yp) CALL Line_To_L12 (x_points = x1_points + tick_length_points, y_points = yp) CALL End_L12_Path (close = .FALSE., stroke = .TRUE., fill = .FALSE.) END DO CALL End_Group() CALL Begin_Group() ! of text labels CALL Set_Fill_or_Pattern (use_pattern = .FALSE., color_name = "foreground") DO j = y_index_low, y_index_high t = j * log10_M_c_tick_interval yp = y1_points + (y2_points - y1_points) * (t - ymin) / (ymax - ymin) IF (log10_M_c_tick_interval == 1.0) THEN ! use integers WRITE (c5, "(I5)") NINT(t) ELSE WRITE (c5, "(F5.2)") t END IF CALL L12_Text (level = 1, x_points = x1_points, y_points = yp, angle_radians = 0.0, & & font_points = 10, lr_fraction = 1.2, ud_fraction = +0.3, & & text = ADJUSTL(c5)) END DO CALL End_Group() CALL L12_Text (level = 1, x_points = x1_points - (text_border_inches * 72), & & y_points = ai_window_yc_points, angle_radians = 1.5708, & & font_points = 11, lr_fraction = 0.5, ud_fraction = +1.0, & & text = "Log10 (Mc / 1 N m)") !Put tick marks and labels along the right vertical axis: corner magnitude minimum_corner_magnitude = (2 * (ymin - 9.05)) / 3 maximum_corner_magnitude = (2 * (ymax - 9.05)) / 3 corner_magnitude_tick_interval = 0.5 y_index_low = Int_Above(minimum_corner_magnitude / corner_magnitude_tick_interval) y_index_high = Int_Below(maximum_corner_magnitude / corner_magnitude_tick_interval) CALL Begin_Group() ! of tick marks CALL Set_Line_Style (width_points = 1.0, dashed = .FALSE.) CALL Set_Stroke_Color("foreground") DO j = y_index_low, y_index_high t = j * corner_magnitude_tick_interval yp = y1_points + (y2_points - y1_points) * (t - minimum_corner_magnitude) / & & (maximum_corner_magnitude - minimum_corner_magnitude) CALL New_L12_Path (level = 1, x_points = x2_points, y_points = yp) CALL Line_To_L12 (x_points = x2_points - tick_length_points, y_points = yp) CALL End_L12_Path (close = .FALSE., stroke = .TRUE., fill = .FALSE.) END DO CALL End_Group() CALL Begin_Group() ! of text labels CALL Set_Fill_or_Pattern (use_pattern = .FALSE., color_name = "foreground") DO j = y_index_low, y_index_high t = j * corner_magnitude_tick_interval yp = y1_points + (y2_points - y1_points) * (t - minimum_corner_magnitude) / & & (maximum_corner_magnitude - minimum_corner_magnitude) IF (corner_magnitude_tick_interval == 1.0) THEN ! use integers WRITE (c5, "(I5)") NINT(t) ELSE WRITE (c5, "(F5.1)") t END IF CALL L12_Text (level = 1, x_points = x2_points, y_points = yp, angle_radians = 0.0, & & font_points = 10, lr_fraction = -0.2, ud_fraction = +0.3, & & text = ADJUSTL(c5)) END DO CALL End_Group() CALL L12_Text (level = 1, x_points = x2_points + (text_border_inches * 72), & & y_points = ai_window_yc_points, angle_radians = 1.5708, & & font_points = 11, lr_fraction = 0.5, ud_fraction = -0.4, & & text = "Corner Magnitude") !Label plot along the top of the contour-map window: CALL Set_Fill_or_Pattern (use_pattern = .FALSE., color_name = "foreground") CALL L12_Text (level = 1, x_points = ai_window_xc_points, & & y_points = y2_points, angle_radians = 0.0, & & font_points = 12, lr_fraction = 0.5, ud_fraction = -0.4, & & text = TRIM(tlab)) !Create inset box for 7 lines of statistical data: CALL Begin_Group() x1_box_points = x1_points + 1.5 * (tick_length_inches * 72) x2_box_points = x1_box_points + (1.4 * 72) y2_box_points = y2_points - 1.5 * (tick_length_inches * 72) y1_box_points = y2_box_points - (1.1 * 72) CALL Set_Line_Style (width_points = 1.0, dashed = .FALSE.) CALL Set_Stroke_Color("foreground") CALL Set_Fill_or_Pattern (use_pattern = .FALSE., color_name = "background") CALL New_L12_Path (level = 1, x_points = x1_box_points, y_points = y1_box_points) CALL Line_To_L12 (x_points = x2_box_points, y_points = y1_box_points) CALL Line_To_L12 (x_points = x2_box_points, y_points = y2_box_points) CALL Line_To_L12 (x_points = x1_box_points, y_points = y2_box_points) CALL Line_To_L12 (x_points = x1_box_points, y_points = y1_box_points) CALL End_L12_Path (close = .TRUE., stroke = .TRUE., fill = .TRUE.) !Insert 7 lines of statistical data into box: likmax(1:6) = 'Lmax =' WRITE(likmax(7:18), '(F12.3)') tabm WRITE (6, 186) likmax WRITE (106, 186) likmax numb(1:11) = 'EQ Number =' WRITE(numb(12:18), '(I7)') nmom WRITE (6, 186) numb WRITE (106, 186) numb gr(1:11) = 'Lmax-GRlik=' WRITE(gr(12:18), '(f7.3)') tabm - tabmgr WRITE (6, 186) gr WRITE (106, 186) gr grt(1:11) = 'Lmax-GRTlk=' WRITE(grt(12:18), '(f7.3)') tabm - tabmgrt WRITE (6, 186) grt WRITE (106, 186) grt magx(1:8) = 'Magmax =' WRITE(magx(9:18), '(F10.3)') amagmax WRITE (6, 186) magx WRITE (106, 186) magx magu(1:8) = 'Mommgx =' WRITE(magu(9:18), '(F10.3)') ammagu WRITE (6, 186) magu WRITE (106, 186) magu betx(1:10) = 'Beta max =' WRITE(betx(11:18), '(F8.3)') betax WRITE (6, 186) betx WRITE (106, 186) betx 186 FORMAT (' ', a18) tx = x1_box_points + (x2_box_points - x1_box_points) * 0.040 tx1 = x1_box_points + (x2_box_points - x1_box_points) * 0.960 ty = y2_box_points - (y2_box_points - y1_box_points) * 0.133 ty1 = y2_box_points - (y2_box_points - y1_box_points) * 0.266 ty2 = y2_box_points - (y2_box_points - y1_box_points) * 0.400 ty3 = y2_box_points - (y2_box_points - y1_box_points) * 0.533 ty4 = y2_box_points - (y2_box_points - y1_box_points) * 0.666 ty5 = y2_box_points - (y2_box_points - y1_box_points) * 0.800 ty6 = y2_box_points - (y2_box_points - y1_box_points) * 0.933 WRITE (6, 166) x1_box_points, x2_box_points, y1_box_points, y2_box_points, tx, ty, ty1 WRITE (106, 166) x1_box_points, x2_box_points, y1_box_points, y2_box_points, tx, ty, ty1 166 FORMAT (' ^^', 8f15.5) CALL Begin_Group() CALL Set_Fill_or_Pattern (use_pattern = .FALSE., color_name = "foreground") !PB CALL Wtstr (tx, ty, likmax, 1, 0, 0) !print max loglikelihood, rick - c CALL L12_Text (level = 1, x_points = tx, y_points = ty, & & angle_radians = 0.0, & & font_points = 10, lr_fraction = 0.0, ud_fraction = 0.0, & & text = likmax(1:6)) CALL L12_Text (level = 1, x_points = tx1, y_points = ty, & & angle_radians = 0.0, & & font_points = 10, lr_fraction = 1.0, ud_fraction = 0.0, & & text = likmax(7:18)) !PB CALL Wtstr (tx, ty1, numb, 1, 0, 0) CALL L12_Text (level = 1, x_points = tx, y_points = ty1, & & angle_radians = 0.0, & & font_points = 10, lr_fraction = 0.0, ud_fraction = 0.0, & & text = numb(1:11)) CALL L12_Text (level = 1, x_points = tx1, y_points = ty1, & & angle_radians = 0.0, & & font_points = 10, lr_fraction = 1.0, ud_fraction = 0.0, & & text = numb(12:18)) !PB CALL Wtstr (tx, ty2, gr, 1, 0, 0) !log ratio tgr / gr 2003 / 1 / 10 CALL L12_Text (level = 1, x_points = tx, y_points = ty2, & & angle_radians = 0.0, & & font_points = 10, lr_fraction = 0.0, ud_fraction = 0.0, & & text = gr(1:11)) CALL L12_Text (level = 1, x_points = tx1, y_points = ty2, & & angle_radians = 0.0, & & font_points = 10, lr_fraction = 1.0, ud_fraction = 0.0, & & text = gr(12:18)) !PB CALL Wtstr (tx, ty3, grt, 1, 0, 0) !log ratio tgr / gr truncated 2003 / CALL L12_Text (level = 1, x_points = tx, y_points = ty3, & & angle_radians = 0.0, & & font_points = 10, lr_fraction = 0.0, ud_fraction = 0.0, & & text = grt(1:11)) CALL L12_Text (level = 1, x_points = tx1, y_points = ty3, & & angle_radians = 0.0, & & font_points = 10, lr_fraction = 1.0, ud_fraction = 0.0, & & text = grt(12:18)) !PB CALL Wtstr (tx, ty4, magx, 1, 0, 0) CALL L12_Text (level = 1, x_points = tx, y_points = ty4, & & angle_radians = 0.0, & & font_points = 10, lr_fraction = 0.0, ud_fraction = 0.0, & & text = magx(1:8)) CALL L12_Text (level = 1, x_points = tx1, y_points = ty4, & & angle_radians = 0.0, & & font_points = 10, lr_fraction = 1.0, ud_fraction = 0.0, & & text = magx(9:18)) !PB CALL Wtstr (tx, ty5, magu, 1, 0, 0) CALL L12_Text (level = 1, x_points = tx, y_points = ty5, & & angle_radians = 0.0, & & font_points = 10, lr_fraction = 0.0, ud_fraction = 0.0, & & text = magu(1:8)) CALL L12_Text (level = 1, x_points = tx1, y_points = ty5, & & angle_radians = 0.0, & & font_points = 10, lr_fraction = 1.0, ud_fraction = 0.0, & & text = magu(9:18)) !PB CALL Wtstr (tx, ty6, betx, 1, 0, 0) CALL L12_Text (level = 1, x_points = tx, y_points = ty6, & & angle_radians = 0.0, & & font_points = 10, lr_fraction = 0.0, ud_fraction = 0.0, & & text = betx(1:10)) CALL L12_Text (level = 1, x_points = tx1, y_points = ty6, & & angle_radians = 0.0, & & font_points = 10, lr_fraction = 1.0, ud_fraction = 0.0, & & text = betx(11:18)) CALL End_Group() CALL End_Group() CALL End_Page() END SUBROUTINE AutoPlot INTEGER FUNCTION Int_Above (x) ! Returns integer equal to, or greater than, x. ! (Note: INT() is different; always truncates toward zero.) IMPLICIT NONE REAL, INTENT(IN) :: x INTEGER :: i REAL :: y i = INT(x) IF (x <= 0.) THEN Int_Above = i ELSE ! x > 0. y = 1.*i IF (y >= x) THEN Int_Above = i ELSE ! most commonly Int_Above = i + 1 END IF END IF END FUNCTION Int_Above INTEGER FUNCTION Int_Below (x) ! Returns integer equal to, or less than, x. ! (Note: INT() is different; always truncates toward zero.) IMPLICIT NONE REAL, INTENT(IN) :: x INTEGER :: i REAL :: y i = INT(x) IF (x >= 0.) THEN Int_Below = i ELSE ! x < 0. y = 1.*i IF (y <= x) THEN Int_Below = i ELSE ! most commonly Int_Below = i - 1 END IF END IF END FUNCTION Int_Below SUBROUTINE Pause() IMPLICIT NONE WRITE (*,"(' Press [Enter]...'\)") READ (*,*) END SUBROUTINE Pause CHARACTER*3 FUNCTION RGB_Munsell(warmth, brightness) !Returns a 3-byte (RGB, or red green blue) pixel color !suitable for packing into a .tiff bitmap array. !Warmth ranges from 0.0 to 1.0, corresponding to the range: ! purples - blues - greens - yellows - browns - reds !Brightness varies from 0.0 to 2.0: ! -at brightness 1.0 there is full color saturation; ! -at brightness below 1.0 some black is mixed in; ! -at brightness above 1.0 some white is mixed in. !Off-scale values of the parameters are treated as equal ! to limiting values. IMPLICIT NONE REAL, INTENT(IN) :: warmth, brightness INTEGER :: blue_integer, green_integer, left_index, red_integer, right_index REAL :: blue_fraction, bright, fraction, green_fraction, & & red_fraction, warm REAL, DIMENSION(3,10) :: Munsell !color (warmth) section: !using Munsell 10-color wheel (copied from Grolier Encyclopedia, ! by eye, so RGB values approximate) DATA Munsell / 255., 80., 255., & ! red-grape & 192., 127., 255., & ! purple & 0., 0., 255., & ! blue & 31., 233., 255., & ! cyan & 0., 215., 103., & ! green & 125., 255., 31., & ! Spring green & 255., 255., 0., & ! yellow & 255., 193., 31., & ! orange & 255., 0., 0., & ! red & 255., 80., 80. / ! lox warm = MAX(0.0, MIN(warmth, 1.0)) left_index = 1 + 9.0 * warm left_index = MIN(9,MAX(1,left_index)) right_index = left_index + 1 fraction = 9.0 * warm - (left_index - 1) fraction = MIN(1.0,MAX(0.0, fraction)) red_fraction = (Munsell(1,left_index) + & & fraction * (Munsell(1,right_index) - Munsell(1,left_index))) / 255.0 green_fraction =(Munsell(2,left_index) + & & fraction * (Munsell(2,right_index) - Munsell(2,left_index))) / 255.0 blue_fraction = (Munsell(3,left_index) + & & fraction * (Munsell(3,right_index) - Munsell(3,left_index))) / 255.0 !brightness section: bright = MAX(0.0, MIN(brightness, 2.0)) IF (bright < 1.0) THEN ! mix with black red_fraction = bright * red_fraction green_fraction = bright * green_fraction blue_fraction = bright * blue_fraction ELSE IF (bright > 1.0) THEN ! mix with white red_fraction = (2.-bright) * red_fraction + bright - 1.0 green_fraction = (2.-bright) * green_fraction + bright - 1.0 blue_fraction = (2.-bright) * blue_fraction + bright - 1.0 ENDIF !expand to range 0-255 and pack into 3 bytes: red_integer = MAX(0,MIN(NINT(255.0* red_fraction),255)) green_integer = MAX(0,MIN(NINT(255.0*green_fraction),255)) blue_integer = MAX(0,MIN(NINT(255.0* blue_fraction),255)) RGB_Munsell = CHAR(red_integer)//CHAR(green_integer)//CHAR(blue_integer) END FUNCTION RGB_Munsell END PROGRAM BetaCorner