SOURCE CODE

This Maximum Entropy program was originally written in BASIC by G.J. Daniell (Department of Physics, Southampton University, UK) and later translated into FORTRAN and adapted for SAS analysis by J.A. Potton. Further modifications have been made by I.D. Culverwell, G.P. Clarke and A.J. Allen (UKAEA Harwell Laboratory,UK) and P.R. Jemian (Northwestern University, USA).

There is only one source code module, MaxSas.For. Compile and link it with the fastest floating point math that you can get your hands on.

Unfortunately, some data storage had to be placed in COMMON because of the limitation of the Language Systems MPW version 1.2.1 FORTRAN compiler for the Apple Macintosh. Because of this compiler’s eccentricacies, there is one compiler-dependent line of code very near the first executable statement. If you use this compiler, un-comment this line so that you get a chance to see the output. (Compiler dependence, ugh!)

As it stands on 7 February 1990, the code will now compile on:

  • Digital Equipment Corporation VAX 11/785,VMS version 5.2
  • Apple Macintosh, Language Systems FORTRAN v. 1.2.1
  • Apple Macintosh, Microsoft (Absoft) FORTRAN v. 2.2 compiler
  • MS-DOS (e.g. IBM-PC), Microsoft FORTRAN v. 5.0

Of course the program RUNS on these computers as well. Quite well!

Most of the comments in the source code have been added by P.R. Jemian. Where they exist, they are usually quite explanatory. Where they do not exist, consult the references of [Skilling, 1984] for the operation of MaxEnt.

Complete listing of MaxSas.for

   1
   2
   3
   4
   5
   6
   7
   8
   9
  10
  11
  12
  13
  14
  15
  16
  17
  18
  19
  20
  21
  22
  23
  24
  25
  26
  27
  28
  29
  30
  31
  32
  33
  34
  35
  36
  37
  38
  39
  40
  41
  42
  43
  44
  45
  46
  47
  48
  49
  50
  51
  52
  53
  54
  55
  56
  57
  58
  59
  60
  61
  62
  63
  64
  65
  66
  67
  68
  69
  70
  71
  72
  73
  74
  75
  76
  77
  78
  79
  80
  81
  82
  83
  84
  85
  86
  87
  88
  89
  90
  91
  92
  93
  94
  95
  96
  97
  98
  99
 100
 101
 102
 103
 104
 105
 106
 107
 108
 109
 110
 111
 112
 113
 114
 115
 116
 117
 118
 119
 120
 121
 122
 123
 124
 125
 126
 127
 128
 129
 130
 131
 132
 133
 134
 135
 136
 137
 138
 139
 140
 141
 142
 143
 144
 145
 146
 147
 148
 149
 150
 151
 152
 153
 154
 155
 156
 157
 158
 159
 160
 161
 162
 163
 164
 165
 166
 167
 168
 169
 170
 171
 172
 173
 174
 175
 176
 177
 178
 179
 180
 181
 182
 183
 184
 185
 186
 187
 188
 189
 190
 191
 192
 193
 194
 195
 196
 197
 198
 199
 200
 201
 202
 203
 204
 205
 206
 207
 208
 209
 210
 211
 212
 213
 214
 215
 216
 217
 218
 219
 220
 221
 222
 223
 224
 225
 226
 227
 228
 229
 230
 231
 232
 233
 234
 235
 236
 237
 238
 239
 240
 241
 242
 243
 244
 245
 246
 247
 248
 249
 250
 251
 252
 253
 254
 255
 256
 257
 258
 259
 260
 261
 262
 263
 264
 265
 266
 267
 268
 269
 270
 271
 272
 273
 274
 275
 276
 277
 278
 279
 280
 281
 282
 283
 284
 285
 286
 287
 288
 289
 290
 291
 292
 293
 294
 295
 296
 297
 298
 299
 300
 301
 302
 303
 304
 305
 306
 307
 308
 309
 310
 311
 312
 313
 314
 315
 316
 317
 318
 319
 320
 321
 322
 323
 324
 325
 326
 327
 328
 329
 330
 331
 332
 333
 334
 335
 336
 337
 338
 339
 340
 341
 342
 343
 344
 345
 346
 347
 348
 349
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363
 364
 365
 366
 367
 368
 369
 370
 371
 372
 373
 374
 375
 376
 377
 378
 379
 380
 381
 382
 383
 384
 385
 386
 387
 388
 389
 390
 391
 392
 393
 394
 395
 396
 397
 398
 399
 400
 401
 402
 403
 404
 405
 406
 407
 408
 409
 410
 411
 412
 413
 414
 415
 416
 417
 418
 419
 420
 421
 422
 423
 424
 425
 426
 427
 428
 429
 430
 431
 432
 433
 434
 435
 436
 437
 438
 439
 440
 441
 442
 443
 444
 445
 446
 447
 448
 449
 450
 451
 452
 453
 454
 455
 456
 457
 458
 459
 460
 461
 462
 463
 464
 465
 466
 467
 468
 469
 470
 471
 472
 473
 474
 475
 476
 477
 478
 479
 480
 481
 482
 483
 484
 485
 486
 487
 488
 489
 490
 491
 492
 493
 494
 495
 496
 497
 498
 499
 500
 501
 502
 503
 504
 505
 506
 507
 508
 509
 510
 511
 512
 513
 514
 515
 516
 517
 518
 519
 520
 521
 522
 523
 524
 525
 526
 527
 528
 529
 530
 531
 532
 533
 534
 535
 536
 537
 538
 539
 540
 541
 542
 543
 544
 545
 546
 547
 548
 549
 550
 551
 552
 553
 554
 555
 556
 557
 558
 559
 560
 561
 562
 563
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585
 586
 587
 588
 589
 590
 591
 592
 593
 594
 595
 596
 597
 598
 599
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609
 610
 611
 612
 613
 614
 615
 616
 617
 618
 619
 620
 621
 622
 623
 624
 625
 626
 627
 628
 629
 630
 631
 632
 633
 634
 635
 636
 637
 638
 639
 640
 641
 642
 643
 644
 645
 646
 647
 648
 649
 650
 651
 652
 653
 654
 655
 656
 657
 658
 659
 660
 661
 662
 663
 664
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697
 698
 699
 700
 701
 702
 703
 704
 705
 706
 707
 708
 709
 710
 711
 712
 713
 714
 715
 716
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741
 742
 743
 744
 745
 746
 747
 748
 749
 750
 751
 752
 753
 754
 755
 756
 757
 758
 759
 760
 761
 762
 763
 764
 765
 766
 767
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
    PROGRAM MaxSAS
    IMPLICIT REAL*8 (A-H,O-Z)
    IMPLICIT INTEGER*4 (I-N)
    CHARACTER*25 ProgVers, EditDate
    PARAMETER ( ProgVers = '3.6 (PRJ)' )
    PARAMETER ( EditDate = '11 February 1992' )
C   Analysis of small-angle scattering data using the technique of
C   entropy maximization.

C   Credits:
C   G.J. Daniell, Dept. of Physics, Southampton University, UK
C   J.A. Potton, UKAEA Harwell Laboratory, UK
C   I.D. Culverwell, UKAEA Harwell Laboratory, UK
C   G.P. Clarke, UKAEA Harwell Laboratory, UK
C   A.J. Allen, UKAEA Harwell Laboratory, UK
C   P.R. Jemian, Northwestern University, USA

C   References:
C   1. J Skilling and RK Bryan; MON NOT R ASTR SOC
C       211 (1984) 111 - 124.
C   2. JA Potton, GJ Daniell, and BD Rainford; Proc. Workshop
C       Neutron Scattering Data Analysis, Rutherford
C       Appleton Laboratory, UK, 1986; ed. MW Johnson,
C       IOP Conference Series 81 (1986) 81 - 86, Institute
C       of Physics, Bristol, UK.
C   3. ID Culverwell and GP Clarke; Ibid. 87 - 96.
C   4. JA Potton, GK Daniell, & BD Rainford,
C       J APPL CRYST 21 (1988) 663 - 668.
C   5. JA Potton, GJ Daniell, & BD Rainford,
C       J APPL CRYST 21 (1988) 891 - 897.

C   This progam was written in BASIC by GJ Daniell and later
C     translated into FORTRAN and adapted for SANS analysis.  It
C     has been further modified by AJ Allen to allow use with a
C     choice of particle form factors for different shapes.  It
C     was then modified by PR Jemian to allow portability between
C     the Digital Equipment Corporation VAX and Apple Macintosh
C     computers.
C   The input data file format is three columns of "Q I dI" which
C     are separated by spaces or tabs.  There is no header line
C     in the input data file.

    PARAMETER (cm2m = 0.01) ! convert cm to m units, but why?
    PARAMETER (MaxPts = 300, MaxBin = 102)
    PARAMETER (isLin = 1, isLog = 2, ioUnit = 1)

C  point-by-point mapping between reciprocal and real space
    COMMON /space1/ grid
    DIMENSION grid(MaxBin,MaxPts)

C  terms used in entropy maximization
    COMMON /space5/ chisq, chtarg, chizer, fSum, blank
    COMMON /space2/ beta, c1, c2, s1, s2
    DIMENSION beta(3), c1(3), c2(3,3), s1(3), s2(3,3)

C  terms used only by subroutine MaxEnt, allocated here to make memory tidy
    COMMON /space3/ ox, z, cgrad, sgrad, xi, eta
    DIMENSION ox(MaxPts), z(MaxPts)
    DIMENSION cgrad(MaxBin), sgrad(MaxBin)
    DIMENSION xi(MaxBin,3), eta(MaxPts,3)

C  space for the plotting frame, allocated here to make memory tidy
C    note the limits: MaxCol <= 100, MaxRow <= 150 (really large screens!)
    PARAMETER (MaxCol = 75, MaxRow = 15)
    PARAMETER (MxC2 = MaxCol+2, MxR2 = MaxRow+2)
    COMMON /space4/ screen, nCol, nRow, nCol2, nRow2
    CHARACTER*1 screen(100, 150)

C  space for main segment arrays
    DIMENSION q(MaxPts), datum(MaxPts), sigma(MaxPts)
    DIMENSION r(MaxBin), f(MaxBin), base(MaxBin), Qty(MaxBin)
    DIMENSION fit(MaxPts), BinWid(MaxBin)
    DIMENSION SkyFit(MaxPts), SkyDis(MaxBin)
    CHARACTER*40 InFile, OutFil
    LOGICAL Yes
    CHARACTER*1 YN, aTab

    DATA one, zero /1.0, 0.0/   ! compiler-independence!
    DATA hrDamp /5.0/   ! model 7&8: sets transition range
    DATA htDamp /0.9/   ! model 7: amount of damping
C  The value "hrDamp" sets the range where the transistion occurs.
C  The value "htDamp" sets the maximum proportion of damping.

C ... Define (initially) the default responses
    DATA iOption    /4/     ! usual form factor for spheres
    DATA Aspect /1.0/       ! particle aspect ratio
    DATA LinLog /isLin/     ! linear binning scale
    DATA n      /40/        ! number of bins
    DATA Dmin, Dmax /8.00, 400.0/   ! particle diameters
    DATA IterMax    /20/        ! maximum number of iterations to try
    DATA RhoSq  /1.0/       ! scattering contrast, x10**28 1/m**4
    DATA fac, err   /1.0, 1.0/  ! scalars for intensity and errors
    DATA qMin, qMax /1.e-8, 100./   ! range to accept
    DATA Bkg    /0.0/       ! intensity to subtract
    DATA sLengt /1.0E-20/   ! rectangular slit-length, 1/A
    DATA SkyBkg /1.0E-6/    ! the so-called "sky background" of [1]

C  Next line for MPW/Language Systems version 1.2.1, Macintosh only
C  Comment this out for other compilers
C  This is the only compiler-dependent line in this source code!!!!!!
C   CALL OutWindowScroll (1000)  ! for 1-line advance screen

    pi = 4. * ATAN(1.)
    aTab = CHAR (9)

C  screen dimension variables for plots, in COMMON /space4/
    nCol = MaxCol
    nRow = MaxRow
    nCol2 = MxC2
    nRow2 = MxR2

    1   WRITE (*,*)
    WRITE (*,*) 'Size distributions from SAS data using the',
     >              ' maximum entropy criterion'
    WRITE (*,*) '       version: ', ProgVers
    WRITE (*,*) '   Last edited: ', EditDate

    CALL GetInf (InFile, OutFil, iOption, Aspect, LinLog,
     >      n, Dmin, Dmax, IterMax, RhoSq, fac, err, qMin,
     >      qMax, Bkg, sLengt, SkyBkg, hrDamp, htDamp)
    IF (InFile .EQ. ' ') STOP

C   Read in the SAS data from the file "InFile"
    WRITE (*,*)  ' Reading from file: ', InFile
        OPEN (UNIT = ioUnit, FILE = InFile, STATUS = 'old')
    DO   j = 1, MaxPts
          READ (ioUnit, *, END = 95) q(j), datum(j), sigma(j)
    END DO
   95   npt=j-1     ! ignore any lines without an explicit EOL mark
        CLOSE (UNIT = ioUnit, STATUS = 'keep') 
    WRITE (*,*) npt, ' points were read from the file'

C   Subtract background, convert to 1/m units and
C   shift for the selected data range
        i = 0
        DO   j = 1, npt
      IF (q(j) .GE. Qmin .AND. q(j) .LE. Qmax) THEN
        i = i + 1
        q(i) = q(j)
        datum(i) = fac * (datum(j)-Bkg) / cm2m
        sigma(i) = fac * err * sigma(j) / cm2m
      END IF
    END DO
    npt = i
    WRITE (*,*) npt, ' points were selected from the data'

C  PRJ: 24 May 1989
C   BinWid: actual radial width of the indexed bin number
C   Step:   radial increment factor (for geometric series)
C   rWid:   radial width (for arithmetic series)
    IF (LinLog .EQ. isLog) THEN ! geometric series
      Step = (Dmax/Dmin)**(1. / FLOAT(n-1)) - 1.
      rWid = 0.
    ELSE                ! arithmetic series
      Step = 0.
      rWid = 0.5*(Dmax - Dmin) / FLOAT(n-1)
    END IF
    r(1) = 0.5 * Dmin
    BinWid(1) = r(1) * Step + rWid
    DO   i = 2, n
      r(i) = r(i-1) + BinWid(i-1)
      BinWid(i) = r(i) * Step + rWid
    END DO

    WRITE (*,*) ' Preparation of the GRID function...'
C  Calculate the form-factor pre-terms
  111   IF (iOption .EQ. 1) THEN    ! Rods, using model of AJ Allen
      alphan1 = 2. * pi * Aspect
      alphan2 = 4. * pi
      preform = alphan1
      sLengt = 0.           ! "pinhole" collimation
    ELSE IF (iOption .EQ. 2) THEN   ! Disks, using model of AJ Allen
      alphan1 = 2. * pi / (Aspect**2)
      alphan2 = 2. * pi
      preform = alphan1
      sLengt = zero
    ELSE IF (iOption .EQ. 3) THEN   ! Globules, using model of AJ Allen
      alphan1 = 4. * pi * Aspect / 3.
      IF (Aspect .LT. 0.99) THEN    ! hamburger-shaped
        sqqt = SQRT (one - Aspect**2)
        argument = (2. - Aspect**2 + 2. * sqqt) / (Aspect**2)
        surchi = (one + Aspect**2 * LOG(argument) / (2.*sqqt) )
     >      / (2. * Aspect)
      ELSE IF (Aspect .GT. 1.01) THEN   ! peanut shaped
        sqqt = SQRT(Aspect**2 - one)
        argument = sqqt / Aspect
        surchi = (one + Aspect**2 * ASIN(argument) / sqqt)
     >      / (2. * Aspect)
      ELSE              ! spheroidal
        surchi = one
      END IF
      alphan2 = 6. * pi * surchi
      preform = alphan1
      sLengt = zero
    ELSE IF (iOption .EQ. 4) THEN   ! Spheres, delta-function
      alphan1 = 4. * pi / 3.
      alphan2 = 6. * pi
      preform = 9. * alphan1
      sLengt = zero
    ELSE IF (iOption .EQ. 5) THEN   ! Spheres, box-distribution
      alphan1 = 4. * pi / 3.    ! This model by PRJ
      alphan2 = 6. * pi
      preform = 48. * pi
      sLengt = zero
    ELSE IF (iOption .EQ. 6) THEN   ! smeared, spheroidal globs
      preform = 4. * Pi / 3.    ! This model by PRJ
      alphan1 = preform
      alphan2 = 6. * Pi
      Cgs = SQRT (3. * Pi)      ! for low-Q region
      Cps = 9. * Pi / 4.        ! for med. high-Q region
      Cp = 9. / 2.          ! for high-Q region
    ELSE IF (iOption .EQ. 7) THEN   ! spheroidal globs, no smearing
      preform = 4. * Pi / 3.    ! This model by PRJ
      alphan1 = preform
      alphan2 = 6. * Pi
      sLengt = zero
    ELSE IF (iOption .EQ. 8) THEN   ! smooth spheres
      preform = 4. * Pi / 3.    ! This model by PRJ
      alphan1 = preform
      alphan2 = 6. * Pi
      sLengt = zero
    END IF

C  alphaN1 is RhoSq * the particle volume
C  alphaN2 is RhoSq * the particle surface area / the particle volume
C   ... and later divided by q**4
    alphan1 = cm2m * alphan1 * rhosq * r(1)**3
    alphan2 = cm2m * alphan2 * rhosq / r(n)
    preform = cm2m * preform * rhosq

    DO   i = 1, n
      rCubed = r(i)**3
      DO   j = 1, npt
        Qr = q(j) * r(i)
        Qr2 = Qr**2
        IF (iOption .EQ. 1) THEN
          QH = q(j) * Aspect * r(i)     ! rod 1/2 - length
          topp = one + 2.*Pi* QH**3 * Qr / (9 * (4 + Qr**2))
     >           + (QH**3 * Qr**4) / 8.
          bott = one + QH**2 * (one + QH**2 * Qr)/9
     >           + (QH**4 * Qr**7) / 16
        ELSE IF (iOption .EQ. 2) THEN
          h = r(i)          ! disk 1/2 - thickness
          Rd = h / Aspect       ! disk radius
          Qh = q(j) * h
          QRd = q(j) * Rd
          topp = one + QRd**3 / (3. + Qh**2)
     >           + (Qh**2 * QRd / 3.)**2
          bott = one + QRd**2 * (one + Qh * QRd**2) / 16
     >           + (Qh**3 * QRd**2 / 3.)**2
        ELSE IF (iOption .EQ. 3) THEN
          topp = one
          bott = one + Qr**2 * (2. + Aspect**2) / 15.
     >           + 2. * Aspect * Qr**4 / (9. * surchi)
        ELSE IF (iOption .EQ. 4) THEN
          topp = (SIN(Qr) - Qr * COS(Qr))**2
          bott = Qr**6
        ELSE IF (iOption .EQ. 5) THEN
          Qj = q(j)
          rP = r(i) + BinWid(i)
          rM = r(i)
          bP = 0.5*rP + (Qj**2)*(rP**3)/6.
     >      + (0.25*(Qj * rP**2) - 0.625/Qj) * SIN (2.*Qj*rP)
     >      + 0.75 * rP * COS (2.*Qj*rP)
          bM = 0.5*rM + (Qj**2)*(rM**3)/6.
     >      + (0.25*(Qj * rM**2) - 0.625/Qj) * SIN (2.*Qj*rM)
     >      + 0.75 * rM * COS (2.*Qj*rM)
          topp = bP - bM
          bott = Qj**6 * (rP**4 - rM**4) * rCubed
        ELSE IF (iOption .EQ. 6) THEN
          rL = r(i) * sLengt
          topp = Cgs
          bott = rL*(one + Qr2/5. + Cgs/Cps * Qr**3)
     >          + Cgs/Cp * Qr**4
        ELSE IF (iOption .EQ. 7) THEN
C  The value "hrDamp" sets the range where the transistion occurs.
C  The value "htDamp" sets the maximum proportion of damping.
C  The weight is a "step" function with a broad edge.
          weight = htDamp * EXP (-Qr2/hrDamp**2) + (one - htDamp)
          topp = 3. * (SIN(Qr) - Qr * COS(Qr)) / Qr**3
          bott = 4.5 / Qr**4    ! bott=<topp**2> for large Qr
          topp = weight * topp**2 + (one-weight) / (one + one/bott)
          bott = one
        ELSE IF (iOption .EQ. 8) THEN  ! like #7 but smoother
          Qr2 = Qr**2
          weight = EXP (-Qr2/hrDamp**2)
          IF (Qr .LE. Pi) THEN
            topp = ((-1./45360.*Qr2+1./840.)*Qr2-1./30.)*Qr2+1./3.
          ELSE
            topp = 0.0
          END IF
          topp = (3*topp)**2
          bott = 4.5 / Qr**4
          topp = weight*topp + (1-weight)/(1 + 1/bott)
          bott = one
        END IF
        grid(i,j) = preform * rCubed * topp / bott
C       factors of 4Pi/3 are already included in "preform"
      END DO
    END DO

C   Attempt to account for scattering from very large and very small 
C   particles by use of the limiting forms of grid(i,j).
    DO   j = 1, npt
      grid(n+1,j) = alphan1 ! next line accounts for a slit-length
      grid(n+2,j) = alphan2 / (q(j)**3 * SQRT(q(j)**2 + sLengt**2))
    END DO

C  Try to solve the problem
C  228  basis = 1.0e-6 / RhoSq      ! Originally was just 1.0e-6
    basis = SkyBkg      ! PRJ, 18.6.90
  228   CALL MaxEnt (n+2,npt, f,datum,sigma, basis,base, max,itermax)

C   "Max" counts the number of iterations inside MAXENT.
C   If Max < IterMax, then the problem has been solved.
    IF (max .GE. itermax) THEN
      WRITE (*,*) ' No convergence! # iter. = ', max
          WRITE (*,*) ' File was: ', InFile
      GO TO 1
    END IF

C   Otherwise, SUCCESS!... so calculate the SAS from the distribution
    CALL opus (n+2, npt, f, fit)
    CALL opus (n+2, npt, base, SkyFit)  ! fit the sky background, too!

C   ... and remove the bin width effect.
C   Also, calculate the total volume fraction, the mode, mean, and
C   standard deviations of the volume and number distributions.
    SumV   = zero
    SumVR  = zero
    SumVR2 = zero
    SumN   = zero
    SumNR  = zero
    SumNR2 = zero
    modeV = 1
    modeN  = 1
    DO   i = 1, n
      size = r(i)
      frac = f(i)
          pVol = 4*Pi/3 * (size * 1.e-8)**3 ! particle volume, cm**3
          IF (iOption .EQ. 1)  pVol = pVol * Aspect ! rods
          IF (iOption .EQ. 2)  pVol = pVol / Aspect ! disks
          IF (iOption .EQ. 3)  pVol = pVol * Aspect ! globs
      amount = (frac - SkyBkg) / pVol       ! number / cm**3
      IF (amount .LT. zero) amount = zero
      f(i) = frac / BinWid(i)
      base(i) = base(i) / BinWid(i)
      Qty(i) = amount / BinWid(i)
      IF (i .GT. 3) THEN            ! ignore 1st few bins
        SumN   = SumN   + amount
        SumNR  = SumNR  + amount * size
        SumNR2 = SumNR2 + amount * size**2
      END IF
      IF (Qty(i) .GT. Qty(modeN)) modeN = i ! get the mode
      SumV   = SumV   + frac
      SumVR  = SumVR  + frac * size
      SumVR2 = SumVR2 + frac * size**2
      IF (f(i) .GT. f(modeV)) modeV = i     ! get the mode
    END DO
    DnMean = 2.0 * SumNR / SumN
    DnSDev = 2.0 * SQRT((SumNR2 / SumN) - (SumNR / SumN)**2)
    DvMean = 2.0 * SumVR / SumV
    DvSDev = 2.0 * SQRT((SumVR2 / SumV) - (SumVR / SumV)**2)

    Entropy = zero
    DO   i = 1, n
      frac = BinWid(i) * f(i) / SumV    ! Skilling & Bryan, eq. 1
      Entropy = Entropy - frac * LOG (frac)
    END DO

C  Show the final distribution, corrected for bin width.

        WRITE (*,*)
    WRITE (*,*) ' Input file: ', InFile
    WRITE (*,*) ' Volume weighted size dist.: V(r)N(r) versus r'
    CALL Plot (n, r, f)

C   Estimate a residual background that remains in the data.
    Sum1 = zero
    Sum2 = zero
    DO  j = 1, npt
      weight = one / (sigma(j)**2)
      Sum1 = Sum1 + weight * (fit(j) -  datum(j))
      Sum2 = Sum2 + weight
    END DO
    shift = Sum1 / Sum2

C  Scale the data back to 1/cm units and calculate Chi-squared
    ChiSq  = zero
    Chi2Bk  = zero
    DO  j = 1, npt
      z(j)  = (datum(j) - fit(j)) / sigma(j) 
      ChiSq   = ChiSq + z(j)**2   
      Chi2Bk  = Chi2Bk + (z(j) + shift/ sigma(j))**2
      datum(j) = cm2m * datum(j)
      sigma(j) = cm2m * sigma(j)
      fit(j) = cm2m * fit(j)
      SkyFit(j) = cm2m * SkyFit(j)
    END DO
    shift = cm2m * shift / fac

    WRITE (*,*) ' standardized residuals vs. point number'
    CALL ResPlt (npt, z)

C  Let the file output begin!

    OPEN (UNIT = ioUnit, FILE=OutFil, STATUS='new')
    WRITE (ioUnit,*) ' Results of maximum entropy analysis of SAS'
    WRITE (ioUnit,*) '     version: ', aTab, ProgVers
    WRITE (ioUnit,*) '      edited: ', aTab, EditDate
    WRITE (ioUnit,*)
    WRITE (ioUnit,*) '  input file: ', aTab, InFile
    WRITE (ioUnit,*) ' output file: ', aTab, OutFil
    WRITE (ioUnit,*)
    WRITE (ioUnit, 35591) 'D, A', aTab, 'f, 1/A',
     >          aTab, 'Bkg f, 1/A', aTab, 'N dD, 1/A/cm^3'
35591   FORMAT (1X, A12, 3(A1, 1X, A15))

    DO  i = 1, n
      WRITE (ioUnit,3559) 2.*r(i), aTab, 0.5*f(i), aTab,
     >                 0.5*Base(i), aTab, 0.5*Qty(i)
    END DO
 3559   FORMAT (1X, F12.2, 3(A1, 1X, 1PE15.5))


    WRITE (ioUnit, 1011) 'Q 1/A', aTab, 'I 1/cm', aTab,
     >          'fit I 1/cm', aTab, 'dI 1/cm', aTab,
     >          'SkyFit 1/cm', aTab, 'z'
 1011   FORMAT (///, A12, 5(1X, A1, A12))

    DO   j = 1, npt
      WRITE (ioUnit,560) q(j), aTab, datum(j), aTab, fit(j),
     >      aTab, sigma(j), aTab, SkyFit(j), aTab, z(j)
    END DO
  560     FORMAT (1PE12.4, 4(A1, E13.5), 1X, A1, 0PF12.6)

    WRITE (ioUnit,3301) aTab, InFile
    WRITE (*,3301) aTab, InFile
 3301   FORMAT (//' Input data: ', A1, A40)

    WRITE (ioUnit,3302) RhoSq
    WRITE (*,3302) RhoSq
 3302   FORMAT (' Contrast = ', F15.7,' x 10^28 m^-4.')

    IF (iOption .EQ. 1) THEN
      WRITE (ioUnit,*) ' rods: dia=D, length=D*', Aspect
      WRITE (*,*) ' rods: dia=D, length=D*', Aspect
    ELSE IF (iOption .EQ. 2) THEN
      WRITE (ioUnit,*) ' disks: thickness=D, dia=D/', Aspect
      WRITE (*,*) ' disks: thickness=D, dia=D/', Aspect
    ELSE IF (iOption .EQ. 3) THEN
      WRITE (ioUnit,*) ' globs: D x D x D*', Aspect
      WRITE (*,*) ' globs: D x D x D*', Aspect
    ELSE IF (iOption .EQ. 4) THEN
      WRITE (ioUnit,*) ' delta-function Spheres: diameter=D'
      WRITE (*,*) ' delta-function Spheres: diameter=D'
    ELSE IF (iOption .EQ. 5) THEN
      WRITE (ioUnit,*) ' box-function Spheres: diameter=D'
      WRITE (*,*) ' box-function Spheres: diameter=D'
    ELSE IF (iOption .EQ. 6) THEN
      WRITE (ioUnit,*) ' slit-smeared spheroidal globs: diameter=D'
      WRITE (*,*) ' slit-smeared spheroidal globs: diameter=D'
      WRITE (ioUnit,*) ' slit-length (1/A) = ', sLengt
      WRITE (*,*) ' slit-length (1/A) = ', sLengt
    ELSE IF (iOption .EQ. 7) THEN
      WRITE (ioUnit,*) ' spheroidal globs: diameter=D'
      WRITE (*,*) ' spheroidal globs: diameter=D'
    ELSE IF (iOption .EQ. 8) THEN
      WRITE (ioUnit,*) ' smooth spheres: diameter=D'
      WRITE (*,*) ' smooth spheres: diameter=D'
    END IF

    WRITE (ioUnit,53303) fac
    WRITE (*,53303) fac
53303   FORMAT (' Data conversion factor to 1/cm = ', 1PE12.5)

    WRITE (ioUnit,63303) err
    WRITE (*,63303) err
63303   FORMAT (' Error scaling factor = ', 1PE12.5)

    IF (LinLog .EQ. isLog) THEN
      WRITE (ioUnit,13304) 'geometric'
      WRITE (*,13304) 'geometric'
    ELSE
      WRITE (ioUnit,13304) 'arithmetic'
      WRITE (*,13304) 'arithmetic'
    END IF
13304   FORMAT (' Histogram bins are distributed in an increasing ',
     >      A10, ' series.')

    WRITE (ioUnit,3304) 'Minimum', Dmin
    WRITE (*,3304) 'Minimum', Dmin
    WRITE (ioUnit,3304) 'Maximum', Dmax
    WRITE (*,3304) 'Maximum', Dmax
 3304   FORMAT (1X, A7, ' particle dimension D = ',F12.2,' A.')

    WRITE (ioUnit,3306) n
    WRITE (*,3306) n
 3306   FORMAT (' Number of histogram bins = ',I4,'.')

    WRITE (ioUnit,3307) itermax
    WRITE (*,3307) itermax
 3307   FORMAT (' Maximum number of iterations allowed = ',I4,'.')

    WRITE (ioUnit,3314) max
    WRITE (*,3314) max
 3314   FORMAT (' Program left MaxEnt routine after ',
     *    I4,' iterations.')

    WRITE (ioUnit,3308) npt
    WRITE (*,3308) npt
 3308   FORMAT (' Target chi-squared (# data points) = ',I5,'.')

    WRITE (ioUnit,3309) ChiSq
    WRITE (*,3309) ChiSq
 3309   FORMAT (' Best value of chi-squared achieved = ',F12.6,'.')

    WRITE (ioUnit, 33091) 'the final', Entropy
    WRITE (*, 33091) 'the final', Entropy
    WRITE (ioUnit, 33091) 'a flat', LOG (FLOAT (n))
    WRITE (*, 33091) 'a flat', LOG (FLOAT (n))
33091   FORMAT (' Entropy of ', A9, ' distribution = ', F12.7,'.')

    WRITE (ioUnit,33101) SumN
    WRITE (*,33101) SumN
33101   FORMAT (' Total particles  = ', 1PE15.5,' per cubic cm.')

    WRITE (ioUnit,3310) SumV
    WRITE (*,3310) SumV
 3310   FORMAT (' Total volume fraction of all scatterers = ',
     *    F15.9,'.')

    WRITE (ioUnit,3311) 'smaller', Dmin, f(n+1)
    WRITE (ioUnit,3311) 'larger',  Dmax, f(n+2)
    WRITE (*,3311) 'smaller', Dmin, f(n+1)
    WRITE (*,3311) 'larger',  Dmax, f(n+2)
 3311   FORMAT (' Volume fraction ',A7,' than ', F12.2,
     *    ' A = ', 1PE13.5,'.')

    WRITE (ioUnit,3411) SkyBkg
    WRITE (*,3411) SkyBkg
 3411   FORMAT (' Sky background (minimum ',
     *    'significant volume fraction) = ', 1PE13.5,'.')

    WRITE (ioUnit,3312) 'Volume', 'mode D value', 2.0 * r(modeV)
    WRITE (*,3312) 'Volume', 'mode D value', 2.0 * r(modeV)
    WRITE (ioUnit,3312) 'Volume', 'mean D value', DvMean
    WRITE (*,3312) 'Volume', 'mean D value', DvMean
    WRITE (ioUnit,3312) 'Volume', 'std. deviation', DvSDev
    WRITE (*,3312) 'Volume', 'std. deviation', DvSDev
    WRITE (ioUnit,3312) 'Number', 'mode D value', 2.0 * r(modeN)
    WRITE (*,3312) 'Number', 'mode D value', 2.0 * r(modeN)
    WRITE (ioUnit,3312) 'Number', 'mean D value', DnMean
    WRITE (*,3312) 'Number', 'mean D value', DnMean
    WRITE (ioUnit,3312) 'Number', 'std. deviation', DnSDev
    WRITE (*,3312) 'Number', 'std. deviation', DnSDev
 3312   FORMAT (1X, A6, '-weighted ', A14, ' = ', F12.5, ' A.')

    WRITE (ioUnit,3313) 'Min', q(1)
    WRITE (*,3313) 'Min', q(1)
    WRITE (ioUnit,3313) 'Max', q(npt)
    WRITE (*,3313) 'Max', q(npt)
 3313   FORMAT (1X, A3,'imum Q-vector = ', 1PE15.7, ' 1/A.')

    WRITE (ioUnit,3315) 'User-specified', Bkg
    WRITE (*,3315) 'User-specified', Bkg
    WRITE (ioUnit,3315) 'Suggested', Bkg - shift
    WRITE (*,3315) 'Suggested', Bkg - shift
 3315   FORMAT (1X, A14, ' background = ', F18.9,' input data units')

    WRITE (ioUnit,*) ' New background should give ChiSq = ', Chi2Bk
    WRITE (*,*) ' New background should give ChiSq = ', Chi2Bk

    CLOSE (UNIT=ioUnit, STATUS='keep')

C  Adjust the background default setting
C  Shift the intensity data just in case the user wants a Stability Check
C  Remember: background shifts down, intensity shifts up
C  Don't forget to put the data back into units of 1/m!
        Bkg = Bkg - shift
    DO   j = 1, npt
      datum(j) = (datum(j) + shift) / cm2m
      sigma(j) = sigma(j) / cm2m
    END DO

    IF (ABS ((Chi2Bk-ChiSq)/FLOAT (npt)) .LE. 0.05) THEN
      WRITE (*,*) ' The change in ChiSquared should be < 5%.'
 4000     WRITE (*, '(X,A,$)') ' Run the Stability Check? (Y/<N>)'
      READ (*,'(A1)') YN
      IF (YN .EQ. 'y'  .OR.  YN .EQ. 'Y') GO TO 228
      IF (YN .NE. ' ' .AND. YN .NE. 'n' .AND. YN .NE. 'N') GO TO 4000
    END IF

    WRITE (*,3200) OutFil
 3200   FORMAT (/,' The program is finished.', /, 
     1    ' The output file is: ', A40)
    GO TO 1

 3199   STOP 
    END


    SUBROUTINE GetInf (InFile, OutFil, iOption, Aspect, LinLog,
     >      nBin, Dmin, Dmax, IterMax, RhoSq, fac, err, qMin,
     >      qMax, Bkg, sLengt, SkyBkg, hrDamp, htDamp)
    IMPLICIT REAL*8 (A-H,O-Z)
    IMPLICIT INTEGER*4 (I-N)
    CHARACTER*40 InFile, OutFil
    PARAMETER (Ro2Max = 1.e6, ItrLim = 200, AbsMax = 1.e3)
    PARAMETER (DiaMin = 1., DiaMax = 1.e6, ErrMax = 1.e6)
    PARAMETER (MaxPts = 300, MaxBin = 102)
    PARAMETER (isLin = 1, isLog = 2)

    1   WRITE (*,'(X,A,$)') ' Input file? <Quit>'
      READ (*, 2) InFile
    2     FORMAT (A40)
      IF (InFile.EQ.' ') RETURN

    3   WRITE (*,'(X,A,$)') ' Output file?'
      READ (*, 2) OutFil
      IF (OutFil .EQ. ' ') GO TO 3
      IF (OutFil .EQ. InFile) GO TO 1

    suggest = qMin
   16   WRITE (*,'(X,A,G,A,$)') ' Minimum q-vector? [1/A] <', 
     >          suggest, '>'
    READ (*, '(F10.0)') qMin
    IF (qMin .LT. 0) GO TO 16
    IF (qMin .EQ. 0) qMin = suggest

    suggest = qMax
   17   WRITE (*,'(X,A,G,A,$)') ' Maximum q-vector? [1/A] <', 
     >          suggest, '>'
    READ (*, '(F10.0)') qMax
    IF (qMax .EQ. 0) qMax = suggest
    IF (qMax .LE. 0) GO TO 17
    IF (qMax .LE. qMin) GO TO 1

    suggest = RhoSq
   13   WRITE (*,'(X,A,G,A,$)') 
     >    ' Scattering contrast? [10^28 m^-4] <', suggest, '>'
    READ (*, '(F10.0)') RhoSq
    IF (RhoSq .EQ. 0) RhoSq = suggest
    IF (RhoSq .LT. 0 .OR. RhoSq .GT. Ro2Max) GO TO 13

    suggest = fac
   14   WRITE (*,'(X,A,G,A,$)') 
     >    ' Factor to convert data to 1/cm? <', suggest, '>'
    READ (*, '(F10.0)') fac
    IF (fac .EQ. 0) fac = suggest
    IF (fac .LE. 0 .OR. fac .GT. AbsMax) GO TO 14

    suggest = err
   15   WRITE (*,'(X,A,G,A,$)') 
     >    ' Error scaling factor? <', suggest, '>'
    READ (*, '(F10.0)') err
    IF (err .EQ. 0) err = suggest
    IF (err .LE. 0 .OR. err .GT. ErrMax) GO TO 15

    suggest = Bkg
   18   WRITE (*,'(X,A,G,A,$)') ' Background? <', suggest, '>'
    READ (*, '(F10.0)') Bkg
    IF (Bkg .EQ. 0) Bkg = suggest

    Last = iOption
    4     WRITE (*,*) '  Select a form model for the scatterer:'
      WRITE (*,*) '  (See the User Guide for complete explanations)'
      WRITE (*,*) ' 1: rods        2: disks       3: globules'
      WRITE (*,*) ' 4: spheres (usual form)       ',
     >          '5: spheres (integrated)'
      WRITE (*,*) ' 6: spheroids (slit-smeared)   ',
     >          '7: spheroidal globs (not smeared)'
      WRITE (*,*) ' 8: smoothed spheres (not smeared)'
      WRITE (*,'(X,A,I3,A,$)') 
     >      '  Which option number?  <', Last, '>'
      READ (*, '(I4)') iOption
      IF (iOption .EQ. 0) iOption = Last
      IF (iOption .LT. 1  .OR.  iOption .GT. 8) GO TO 4

    suggest = Aspect
    6   IF (iOption .GE. 1  .AND.  iOption .LE. 3) THEN
      WRITE (*,*) ' AR = Aspect Ratio, useful ranges are indicated'
      IF (iOption .EQ. 1) THEN
        WRITE (*,*) ' diameter D, length D * AR, AR > 5'
      ELSE IF (iOption .EQ. 2) THEN
        WRITE (*,*) ' thickness D, diameter D / AR, AR < 0.2'
      ELSE IF (iOption .EQ. 3) THEN
        WRITE (*,*) ' D x D x D * AR, 0.3 < AR < 3'
      END IF
      WRITE (*,'(X,A,G,A,$)') 
     >      ' Aspect ratio?  <', suggest, '>'
      READ (*,'(F10.0)') Aspect
      IF (Aspect .EQ. 0) Aspect = suggest
      IF (Aspect .LT. 0) GO TO 6
    END IF

    suggest = sLengt
   61   IF (iOption .EQ. 6) THEN
      WRITE (*,'(X,A,G,A,$)') 
     >      ' Slit-smeared globs.  Slit-length [1/A]? <', 
     >      suggest, '>'
      READ (*,'(F10.0)') sLengt
      IF (sLengt .EQ. 0) sLengt = suggest
      IF (sLengt .LT. 0) GO TO 61
    END IF

    suggest = htDamp
   62   IF (iOption .EQ. 7) THEN
      WRITE (*,'(X,A,G,A,$)') 
     >      ' spheroidal globs.  fraction of standard function? <', 
     >      suggest, '>'
      READ (*,'(F10.0)') htDamp
      IF (htDamp .EQ. 0) htDamp = suggest
      IF (htDamp .LT. 0) GO TO 62
      IF (htDamp .GT. 1) GO TO 62
    END IF

    suggest = hrDamp
   63   IF (iOption .EQ. 7  .OR.  iOption .EQ. 8) THEN
      WRITE (*,'(X,A,G,A,$)') 
     >      ' smoothed spheres.  Onset Qr value? <', 
     >      suggest, '>'
      READ (*,'(F10.0)') hrDamp
      IF (hrDamp .EQ. 0) hrDamp = suggest
      IF (hrDamp .LT. 0) GO TO 63
    END IF

    Last = LinLog
    7     WRITE (*,'(X,A,I2,A,$)') 
     >      ' Bin step scale? (1=Linear, 2=Log) <', Last, '>'
    READ (*, '(I4)') LinLog
    IF (LinLog .EQ. 0) LinLog = Last
    IF (LinLog .NE. isLin  .AND. LinLog .NE. isLog) GO TO 7

    Last = nBin
    8     WRITE (*,'(X,A,I4,A,$)') 
     >      ' Number of histogram bins? <', Last, '>'
    READ (*, '(I4)') nBin
    IF (nBin .EQ. 0) nBin = Last
    IF (nBin .LT. 2 .OR. nBin .GT. (MaxBin-2)) GO TO 8

    suggest = Dmax
    9     WRITE (*,'(X,A,G,A,$)') 
     >      ' Maximum value of D? [A] <', suggest, '>'
    READ (*, '(F10.0)') Dmax
    IF (Dmax .EQ. 0) Dmax = suggest
    IF (Dmax .LT. nBin*DiaMin .OR. Dmax .GE. DiaMax) GO TO 9

    Suggest = Dmax / FLOAT (nBin)
   11     WRITE (*,'(X,A,G,A,$)') 
     >      ' Minimum value of D? [A] <', suggest, '>'
    READ (*, '(F10.0)') Dmin
    IF (Dmin .EQ. 0) Dmin = suggest
    IF (Dmin .GE. DMax .OR. Dmin .LT. DiaMin) GO TO 1

    IF (IterMax .GT. ItrLim) IterMax = ItrLim
    Last = IterMax
   12     WRITE (*,'(X,A,I4,A,$)') 
     >      ' Maximum number of iterations? <', Last, '>'
    READ (*, '(I4)') IterMax
    IF (IterMax .EQ. 0) IterMax = Last
    IF (IterMax .LT. 0 .OR. IterMax .GT. ItrLim) GO TO 12

    Suggest = SkyBkg
   21     WRITE (*,'(X,A,G,A,$)') 
     >      ' Sky background? (positive) <', Suggest, '>'
    READ (*, '(F10.0)') SkyBkg
    IF (SkyBkg .LT. 0) GO TO 21
    IF (SkyBkg .EQ. 0) SkyBkg = Suggest ! keep default

    RETURN
    END


    SUBROUTINE opus(n,npt,x,ox) ! solution-space -> data-space
    IMPLICIT REAL*8 (A-H,O-Z)
    IMPLICIT INTEGER*4 (I-N)
    PARAMETER (MaxPts=300, MaxBin=102)
    COMMON /space1/ grid
    DIMENSION x(MaxBin), grid(MaxBin,MaxPts), ox(MaxPts)
    DO   j = 1, npt
      sum = 0.
      DO   i = 1, n
       sum = sum + x(i) * grid(i,j)
      END DO
      ox(j) = sum
    END DO
    RETURN
    END


    SUBROUTINE tropus(n,npt,ox,x)   ! data-space -> solution-space
    IMPLICIT REAL*8 (A-H,O-Z)
    IMPLICIT INTEGER*4 (I-N)
    PARAMETER (MaxPts=300, MaxBin=102)
    COMMON /space1/ grid
    DIMENSION x(MaxBin), grid(MaxBin,MaxPts), ox(MaxPts)
    DO   i = 1, n
      sum = 0.
      DO   j = 1, npt
        sum = sum + ox(j) * grid(i,j)
      END DO
      x(i) = sum
    END DO
    RETURN
    END


    SUBROUTINE MaxEnt(n,npt, f,datum,sigma, flat,base,iter,itermax)
    IMPLICIT REAL*8 (A-H,O-Z)
    IMPLICIT INTEGER*4 (I-N)
    PARAMETER (MaxPts=300, MaxBin=102)
    DIMENSION f(MaxBin), datum(MaxPts), sigma(MaxPts)
    DIMENSION base(MaxBin)

    COMMON /space1/ grid
    DIMENSION grid(MaxBin,MaxPts)

    COMMON /space5/ chisq, chtarg, chizer, fSum, blank
    COMMON /space2/ beta, c1, c2, s1, s2
    PARAMETER (m = 3)       ! number of search directions
    DIMENSION beta(m), c1(m), c2(m,m), s1(m), s2(m,m)

    COMMON /space3/ ox, z, cgrad, sgrad, xi, eta
    DIMENSION ox(MaxPts), z(MaxPts)
    DIMENSION cgrad(MaxBin), sgrad(MaxBin)
    DIMENSION xi(MaxBin,3), eta(MaxPts,3)

    PARAMETER (TstLim = 0.05)   ! for convergence
    DATA one, zero /1.0, 0.0/   ! compiler-independence!

        WRITE (*,*) ' MaxEnt routine beginning ...'

    chizer = FLOAT(npt)
    chtarg = chizer
    blank = flat
    exp1 = EXP(one)

    IF (blank .EQ. zero) THEN
      DO  i = 1, n
        blank = blank + base(i)
        f(i) = base(i)  ! given initial distribution
      END DO
      blank = blank / FLOAT(n)
      WRITE (*,*) ' Average of BASE = ', blank
    ELSE
      WRITE (*,*) ' Setting BASE constant at ', blank
      DO  i = 1, n
        base(i) = blank
        f(i) = blank    ! featureless initial distribution
      END DO
    ENDIF

    iter = 0
    6   iter = iter + 1     ! The iteration loop begins here!
    CALL opus (n, npt, f, ox)   ! calc. the model intensity from "f"
    chisq = zero
    DO  j = 1, npt
      a = (ox(j) - datum(j)) / sigma(j)
      chisq = chisq + a**2
      ox(j) = 2. * a / sigma(j)
    END DO
    CALL tropus(n,npt,ox,cgrad) ! cGradient = Grid * ox
    test = zero ! mismatch between entropy and ChiSquared gradients
    snorm = zero    ! entropy term
    cnorm = zero    ! ChiSqr term
    tnorm = zero    ! norm for the gradient term TEST
    fSum = zero ! find the sum of the f-vector
    DO   i = 1, n
      fSum = fSum + f(i)
      sgrad(i) = -LOG(f(i)/base(i)) / (base(i)*exp1)
      snorm = snorm + f(i) * sgrad(i)**2
      cnorm = cnorm + f(i) * cgrad(i)**2
      tnorm = tnorm + f(i) * sgrad(i) * cgrad(i)
    END DO
    snorm = SQRT(snorm)
    cnorm = SQRT(cnorm)
    a = one
    b = one / cnorm
    IF (iter .GT. 1) THEN
      test = SQRT(0.5*(one-tnorm/(snorm*cnorm)))
      a = 0.5 / (snorm * test)
      b = 0.5 * b / test
    ENDIF
    DO  i = 1, n
      xi(i,1) = f(i) * cgrad(i) / cnorm
      xi(i,2) = f(i) * (a * sgrad(i) - b * cgrad(i))
    END DO
    CALL opus (n,npt,xi(1,1),eta(1,1))
    CALL opus (n,npt,xi(1,2),eta(1,2))
    DO  j = 1, npt
      ox(j) = eta(j,2) / (sigma(j)**2)
    END DO
    CALL tropus (n,npt,ox,xi(1,3))
    a = zero
    DO  i = 1, n
      b = f(i) * xi(i,3)
      a = a + b * xi(i,3)
      xi(i,3) = b
    END DO
    a = one / SQRT(a)
    DO  i = 1, n
      xi(i,3) = a * xi(i,3)
    END DO
    CALL opus (n,npt,xi(1,3),eta(1,3))
    DO  k = 1, m
      s1(k) = zero
      c1(k) = zero
      DO  i = 1, n
        s1(k) = s1(k) + xi(i,k) * sgrad(i)
        c1(k) = c1(k) + xi(i,k) * cgrad(i)
      END DO
      c1(k) = c1(k) / chisq
    END DO
    DO  k = 1, m
      DO  l = 1, k
        s2(k,l) = zero
        c2(k,l) = zero
        DO  i = 1, n
          s2(k,l) = s2(k,l) - xi(i,k) * xi(i,l) / f(i)
        END DO
        DO  j = 1, npt
          c2(k,l) = c2(k,l) + eta(j,k) * eta(j,l) / (sigma(j)**2)
        END DO
        s2(k,l) = s2(k,l) / blank
        c2(k,l) = 2. * c2(k,l) / chisq
      END DO
    END DO
    c2(1,2) = c2(2,1)
    c2(1,3) = c2(3,1)
    c2(2,3) = c2(3,2)
    s2(1,2) = s2(2,1)
    s2(1,3) = s2(3,1)
    s2(2,3) = s2(3,2)
    beta(1) = -0.5 * c1(1) / c2(1,1)
    beta(2) = zero
    beta(3) = zero
    IF (iter .GT. 1) CALL Move(m)

C  Modify the current distribution (f-vector)
    fSum = zero     ! find the sum of the f-vector
    fChange = zero      ! and how much did it change?
    DO  i = 1, n
      df = beta(1)*xi(i,1)+beta(2)*xi(i,2)+beta(3)*xi(i,3)
      IF (df .LT. -f(i)) df = 0.001 * base(i) - f(i)    ! a patch
      f(i) = f(i) + df      ! adjust the f-vector
      fSum = fSum + f(i)
      fChange = fChange + df
    END DO

    s = zero
    DO   i = 1, n
      temp = f(i) / fSum        ! fraction of f(i) in this bin
      s = s - temp * LOG (temp) ! from Skilling and Bryan, eq. 1
    END DO

        CALL opus (n, nPt, f, z)    ! model the data-space from f(*)
    ChiSq = zero            ! get the new ChiSquared
        DO   j = 1, nPt
          z(j) = (datum(j) - z(j)) / sigma(j)   ! the residuals
      ChiSq = ChiSq + z(j)**2   ! report this ChiSq, not the one above
    END DO

  300   IF ( MOD(iter, 5) .EQ. 0 ) THEN
      WRITE (*,*)
      WRITE (*,*) ' Residuals'
      CALL ResPlt (npt, z)

      WRITE (*,*)
      WRITE (*,*) ' Distribution'
      CALL BasPlt (n, f, base)
    END IF

    WRITE (*,*) ' #', iter, ' of ', itermax, ',  n  = ', npt
    WRITE (*,200) test, s
    WRITE (*,201) 'target',SQRT(chtarg/npt), 'now',SQRT(chisq/npt)
    WRITE (*,202) 'sum', fSum, ' % change', 100.*fChange/fSum
  200   FORMAT (' test = ', F9.5, ',  Entropy = ', F12.7)
  201   FORMAT (' SQRT((Chi^2)/n):', A8,' = ', F12.8,A10,' = ', F12.8)
  202   FORMAT ('        f-vector:', A8,' = ', F12.8,A10,' = ', F12.8)

C  See if we have finished our task.
    IF (ABS(chisq/chizer-one) .LT. 0.01) THEN  ! hardest test first
      IF (test .LT. TstLim) THEN        ! same solution gradient?
C       We've solved it but now must check for a bizarre condition.
C       Calling routine says we failed if "iter = iterMax".
C       Let's increment iterMax so (maybe) this doesn't happen.
        IF (iter .EQ. iterMax) iterMax = iterMax + 1
        RETURN
      END IF
    END IF
    IF (iter .LT. iterMax) GO TO 6

C  Ask for more time to finish the job.
    WRITE (*,*)
    WRITE (*,*) ' Maximum iterations have been reached.'
 2001   WRITE (*,*) ' How many more iterations? <none>'
    READ (*,'(I4)') more
    IF (more .LT. 0) GO TO 2001
    IF (more .EQ. 0) RETURN
    iterMax = iterMax + more
    GO TO 6
    END


    SUBROUTINE Move(m)
    IMPLICIT REAL*8 (A-H,O-Z)
    IMPLICIT INTEGER*4 (I-N)
    PARAMETER ( MxLoop = 500 )  ! for no solution
    PARAMETER ( Passes = 1.e-3 )    ! convergence test
    COMMON /space5/ chisq, chtarg, chizer, fSum, blank
    COMMON /space2/ beta, c1, c2, s1, s2
    DIMENSION beta(3), c1(3), c2(3,3), s1(3), s2(3,3)
    DATA one, zero /1.0, 0.0/   ! compiler-independence!
    a1 = zero           ! lower bracket  "a"
    a2 = one            ! upper bracket of "a"
        cmin = ChiNow (a1, m)
    IF (cmin*chisq .GT. chizer) ctarg = 0.5*(one + cmin)
    IF (cmin*chisq .LE. chizer) ctarg = chizer/chisq
    f1 = cmin - ctarg
    f2 = ChiNow (a2,m) - ctarg
    DO   loop = 1, MxLoop
      anew = 0.5 * (a1+a2)      ! choose a new "a"
      fx = ChiNow (anew,m) - ctarg
      IF (f1*fx .GT. zero) a1 = anew
      IF (f1*fx .GT. zero) f1 = fx
      IF (f2*fx .GT. zero) a2 = anew
      IF (f2*fx .GT. zero) f2 = fx
      IF (abs(fx) .LT. Passes) GO TO 2
    END DO

C  If the preceding loop finishes, then we do not seem to be converging.
C   Stop gracefully because not every computer uses control-C (etc.)
C   as an exit procedure.
    WRITE (*,*) ' Loop counter = ', MxLoop
    PAUSE ' No convergence in alpha chop (MOVE).  Press return ...'
    STOP ' Program cannot continue.'

    2   w = Dist (m)
    IF (w .GT. 0.1*fSum/blank) THEN
      DO  k = 1, m
        beta(k) = beta(k) * SQRT(0.1 * fSum/(blank * w))
      END DO
    END IF
    chtarg = ctarg * chisq
    RETURN
    END


    REAL*8 FUNCTION Dist (m)
    IMPLICIT REAL*8 (A-H,O-Z)
    IMPLICIT INTEGER*4 (I-N)
    COMMON /space5/ chisq, chtarg, chizer, fSum, blank
    COMMON /space2/ beta, c1, c2, s1, s2
    DIMENSION beta(3), c1(3), c2(3,3), s1(3), s2(3,3)
    DATA one, zero /1.0, 0.0/   ! compiler-independence!
    w = zero
    DO   k = 1, m
      z = zero
      DO   l = 1, m
        z = z - s2(k,l) * beta(l)
      END DO
      w = w + beta(k) * z
    END DO
    Dist = w
    RETURN
    END


    REAL*8 FUNCTION ChiNow(ax,m)
    IMPLICIT REAL*8 (A-H,O-Z)
    IMPLICIT INTEGER*4 (I-N)
    COMMON /space5/ chisq, chtarg, chizer, fSum, blank
    COMMON /space2/ beta, c1, c2, s1, s2
    DIMENSION beta(3), c1(3), c2(3,3), s1(3), s2(3,3)
    DIMENSION a(3,3), b(3)
    DATA one, zero /1.0, 0.0/   ! compiler-independence!
    bx = one - ax
    DO   k = 1, m
      DO   l = 1, m
        a(k,l) = bx * c2(k,l)  -  ax * s2(k,l)
      END DO
      b(k) = -(bx * c1(k)  -  ax * s1(k))
    END DO
    CALL ChoSol(a,b,m,beta)
    w = zero
    DO   k = 1, m
      z = zero
      DO   l = 1, m
        z = z + c2(k,l) * beta(l)
      END DO
      w = w + beta(k) * (c1(k) + 0.5 * z)
    END DO
    ChiNow = one +  w
    RETURN
    END


    SUBROUTINE ChoSol(a, b, n, beta)
    IMPLICIT REAL*8 (A-H,O-Z)
    IMPLICIT INTEGER*4 (I-N)
    DIMENSION fl(3,3), a(3,3), bl(3), b(3), beta(3)
    DATA one, zero /1.0, 0.0/   ! compiler-independence!
    IF (a(1,1) .LE. zero) THEN
      WRITE (*,*) ' Fatal error in CHOSOL: a(1,1) = ', a(1,1)
      PAUSE ' Press <RETURN> to end program ...'
      STOP ' Program cannot continue.'
    END IF
    fl(1,1) = SQRT(a(1,1))
    DO   i = 2, n
      fl(i,1) = a(i,1) / fl(1,1)
      DO   j = 2, i
        z = zero
        DO   k = 1, j-1
          z = z + fl(i,k) * fl(j,k)
        END DO
        z = a(i,j) - z
        IF (j .EQ. i) fl(i,j) = SQRT(z)
        IF (j .NE. i) fl(i,j) = z / fl(j,j)
      END DO
    END DO
    bl(1) = b(1) / fl(1,1)
    DO   i=2, n
      z = zero
      DO   k = 1, i-1
        z = z + fl(i,k) * bl(k)
      END DO
      bl(i) = (b(i) - z) / fl(i,i)
    END DO
    beta(n) = bl(n) / fl(n,n)
    DO   i1 = 1, n-1
      i = n - i1
      z = zero
      DO   k = i+1, n
        z = z + fl(k,i) * beta(k)
      END DO
      beta(i) = (bl(i) - z) / fl(i,i)
    END DO
    RETURN
    END


    SUBROUTINE ResPlt (n, x)
C   Draw a plot of the standardized residuals on the screen.
C   Mark the rows of + and - one standard deviation.
    IMPLICIT REAL*8 (A-H,O-Z)
    IMPLICIT INTEGER*4 (I-N)
    DIMENSION x(1)
    CHARACTER*1 Blank, Symbol, hBordr, vBordr, resSym
    PARAMETER (Blank = ' ', Symbol = 'O', resSym = '=')
    PARAMETER (hBordr = '-', vBordr = '|')

    COMMON /space4/ screen, MaxCol, MaxRow, MxC2, MxR2
    CHARACTER*1 screen(100, 150)

        IF (n .LT. 2) RETURN    ! not enough data

C  Find out how many points to pack per column and how many columns
    nPack = 1 + INT(FLOAT (n) / MaxCol - 1./n)
    nCol = INT((n - 1./n)/nPack + 1)

C  prepare the "screen" for drawing
    DO   j = 1, nCol + 2
      DO   i = 1, MxR2
        screen(i,j) = Blank
      END DO
    END DO
    DO   i = 2, nCol + 1
      screen(MxR2,i) = hBordr
      screen(1,i) = hBordr
    END DO
    DO   i = 2, MaxRow + 1
      screen(i,nCol+2) = vBordr
      screen(i,1) = vBordr
    END DO

C  get the data limits
        xMax = 1.
    xMin = -1.
        DO  i = 1, n
      IF (x(i) .GT. xMax) xMax = x(i)
      IF (x(i) .LT. xMin) xMin = x(i)
    END DO
    RowDel = (MaxRow - 1) / (xMax - xMin)

C  show the standard deviation bars
    mPlus = 1 + INT((1 - xMin)*RowDel + 1)
    mMinus = 1 + INT((-1 - xMin)*RowDel + 1)
    DO   i = 2, nCol + 1
      screen(mMinus,i) = resSym
      screen(mPlus,i) = resSym
    END DO

C  draw the data (overdrawing the residuals bars if necessary)
    DO   i = 1, n
      mCol = 1 + INT((i - 1./n)/nPack + 1)      ! addressing function
      mRow = 1 + INT((x(i) - xMin)*RowDel + 1)  ! +1 for the plot frame
      screen(mRow, mCol) = Symbol
    END DO

C  convey the "screen" to the default output
    WRITE (*,*) nPack, ' point(s) per column'
    WRITE (*,*) 1./RowDel, ' standard deviations per row'
    DO   i = MxR2, 1, -1
      WRITE (*,*) (screen(i,j), j = 1, nCol + 2)
    END DO

    RETURN
    END


    SUBROUTINE BasPlt (n, x, basis)
C   Draw a plot of some data with reference to a basis line on the plot.
C   The basis is that line below which the data is not meaningful.
    IMPLICIT REAL*8 (A-H,O-Z)
    IMPLICIT INTEGER*4 (I-N)
    DIMENSION x(1), basis(1)
    CHARACTER*1 Blank, Symbol, hBordr, vBordr, BasSym
    PARAMETER (Blank = ' ', Symbol = 'O', BasSym = '=')
    PARAMETER (hBordr = '-', vBordr = '|')

    COMMON /space4/ screen, MaxCol, MaxRow, MxC2, MxR2
    CHARACTER*1 screen(100, 150)

        IF (n .LT. 2) RETURN    ! not enough data

C  Find out how many points to pack per column and how many columns
    nPack = 1 + INT(FLOAT (n) / MaxCol - 1./n)
    nCol = INT((n - 1./n)/nPack + 1)

C  prepare the "screen" for drawing
    DO  j = 1, nCol + 2
      DO  i = 1, MxR2
        screen(i,j) = Blank
      END DO
    END DO
    DO  i = 2, nCol + 1
      screen(MxR2,i) = hBordr
      screen(1,i) = hBordr
    END DO
    DO  i = 2, MaxRow + 1
      screen(i,nCol+2) = vBordr
      screen(i,1) = vBordr
    END DO

C  get the data limits
        xMax = x(1)
    xMin = xMax
        DO  i = 1, n
      IF (x(i) .GT. xMax) xMax = x(i)
      IF (x(i) .LT. xMin) xMin = x(i)
      IF (basis(i) .GT. xMax) xMax = basis(i)
      IF (basis(i) .LT. xMin) xMin = basis(i)
    END DO
    RowDel = (MaxRow - 1) / (xMax - xMin)


C  draw the data (overdrawing the basis bars if necessary)
    DO  i = 1, n
      mCol = 1 + INT((i - 1./n)/nPack + 1)      ! addressing function
      mRow = 1 + INT((basis(i) - xMin)*RowDel + 1)  ! basis
      screen(mRow, mCol) = basSym
      mRow = 1 + INT((x(i) - xMin)*RowDel + 1)  ! data
      screen(mRow, mCol) = Symbol
    END DO

C  convey the "screen" to the default output
    WRITE (*,*) nPack, ' point(s) per column'
    WRITE (*,*) 1./RowDel, ' units per row'
    DO  i = MxR2, 1, -1
      WRITE (*,*) (screen(i,j), j = 1, nCol + 2)
    END DO

    RETURN
    END


        SUBROUTINE Plot (n,x,y)
C   Make a scatter plot on the default display device (UNIT=*).
C   MaxRow and MaxCol correspond to the display dimensions.
    IMPLICIT REAL*8 (A-H,O-Z)
    IMPLICIT INTEGER*4 (I-N)
    DIMENSION x(1), y(1)
    CHARACTER*1 Blank, Symbol, hBordr, vBordr
    PARAMETER (Blank = ' ', Symbol = 'O')
    PARAMETER (hBordr = '-', vBordr = '|')

    COMMON /space4/ screen, MaxCol, MaxRow, MxC2, MxR2
    CHARACTER*1 screen(100, 150)

        IF (n .LT. 2) RETURN    ! not enough data

C  prepare the "screen" for drawing
    DO   j = 1, MxC2
      DO   i = 1, MxR2
        screen(i,j) = Blank
      END DO
    END DO
    DO   i = 2, MaxCol+1
      screen(MxR2,i) = hBordr
      screen(1,i) = hBordr
    END DO
    DO   i = 2, MaxRow+1
      screen(i,MxC2) = vBordr
      screen(i,1)  = vBordr
    END DO

C  get the data limits
    xMin = x(1)
        xMax = x(1)
    yMin = y(1)
        yMax = y(1)
        DO  i = 2, n
      IF (x(i).GT.xMax) xMax=x(i)
      IF (x(i).LT.xMin) xMin=x(i)
      IF (y(i).GT.yMax) yMax=y(i)
      IF (y(i).LT.yMin) yMin=y(i)
    END DO
    ColDel = (MaxCol - 1) / (xMax - xMin)
    RowDel = (MaxRow - 1) / (yMax - yMin)

C  data scaling functions are offset by +1 for plot frame
    DO   i = 1, n
      mCol = 1 + INT((x(i) - xMin)*ColDel + 1)
      mRow = 1 + INT((y(i) - yMin)*RowDel + 1)
      screen(mRow, mCol) = Symbol
    END DO

C  convey the "screen" to the default output
    WRITE (*,*) 1./ColDel, ' units per column'
    WRITE (*,*) 1./RowDel, ' units per row'
    DO   i = MaxRow + 2, 1, -1
      WRITE (*,*) (screen(i,j), j = 1, MaxCol + 2)
    END DO
        RETURN
        END

Table Of Contents

Previous topic

AN ALTERNATE TEST DISTRIBUTION: REVERSE.SAS (by P.R. Jemian)

This Page