Linux-Unix program

Development Platform:

Fortran

contin.for：Code Content

C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0003
C CONTIN. MAIN SUBPROGRAM. 0004
c
C Download the Users Manual and other essential documentation from
c http://S-provencher.COM
c
c 06-Dec-01: 500 light scattering form factors also allowed (increased from the
c unnecessarily small 50).
C
C 20-Jan-99: DIMENSION specifications increased to handle up to 8192 data
C points and about 500 grid points. Section 3.5 of the manual
C explains how you can change these, if memory is limited.
C If computation time is critical and memory limited, then
C reducing the maximum number of grid points as much as
C possible might help a little.
C Increasing the DIMENSION specifications for more than 8192
C data points is no problem.
C More than 500 grid points is not recommended, because of
C possible problems with numerical stability.
C With difficult inversions, such as Laplace transforms, NG,
C the number of grid points, should probably not exceed 200.
C However, NG is only an input parameter; you don't have to
C reduce the DIMENSION specifications or make any changes to
C the code.
C
C In each run, all the DIMENSION specifications are tested. Thus,
C you can simply try running, and see if CONTIN aborts with a
C diagnostic telling you to increase a DIMENSION specification.
C If it does not abort, then everything is OK.
C
C==============================================================================
C Hints for inverting noisy Laplace transforms:
C =============================================
C With Laplace inversions, you might use the following Control
C Parameters:
C NONNEG=T (if the solution is known to be nonnegative; without
C this constraint the resolution is extremely poor.)
C IUSER(10)=2
C GMNMX(1) about 0.1/(maximum time value in your data of
C signal vs. time)
C GMNMX(2) about 4/(time-spacing between data points at the
C shortest times)
C NG about 12*log10[GMNMX(2)/GMNMX(1)]
C IFORMY, NINTT, etc., as appropriate for your data.
C
C The CHOSEN SOLUTION gives you a conservative (smooth) estimate
C of a possible continuous distribution of exponentials.
C
C The Reference Solution (the solution with the smallest ALPHA)
C gives you the optimal analysis as a discrete sum of
C exponentials. The number of discrete exponentials is the
C number of peaks. Each amplitude is given by MOMEMT(0) for
C the corresponding peak. The decay rate constant is given by
C MOMENT(1)/MOMENT(0).
C
C Choose the solution with the smallest ALPHA that has the same
C number of peaks as the CHOSEN SOLUTION. This solution has
C less smoothing bias (smaller ALPHA), but still has about
C the same complexity (number of peaks) as the CHOSEN SOLUTION.
C===============================================================================
C
C 5-Mar-94: COMMON blocks /MS1/-/MS6/ were added to simplify compilation on
C PCs that have restricted easily useable memory.
C
C With MS-DOS, you will probably have to break this file into 3 or 4
C segments for compiling and then LINK these object files together.
C You will probably also have to drastically reduce the DIMENSION
C specifications for the maximum number of grid points (down to
C about 50); Section 3.5 of the manual tells you how to do this.
C
C
C FOR THE REGULARIZED SOLUTION OF LINEAR ALGEBRAIC AND 0005
C LINEAR FREDHOLM INTEGRAL EQUATIONS OF THE FIRST KIND, WITH 0006
C OPTIONS FOR PEAK CONSTRAINTS AND LINEAR EQUALITY AND INEQUALITY 0007
C CONSTRAINTS. 0008
C REFERENCES - S.W. PROVENCHER (1982) COMPUT. PHYS. COMMUN., VOL. 27, 0009
C PAGES 213-227, 229-242. 0010
C (1984) CONTIN USERS MANUAL (EMBL 0011
C TECHNICAL REPORT DA07). 0012
C
C As a new user you should notify S.W. Provencher at:
C sp@S-provencher.COM
C 0027
C----------------------------------------------------------------------- 0028
C CALLS SUBPROGRAMS - BLOCK DATA, INIT, INPUT, SETGRD, USERSI, 0029
C WRITYT, USERSX, USERNQ, SETNNG, ANALYZ, SETWT 0030
C WHICH IN TURN CALL - STORIN, READYT, ERRMES, WRITIN, USERIN, 0031
C USERGR, CQTRAP, USERTR, USEREX, RGAUSS, RANDOM, SETSCA, SEQACC, 0032
C H12, GETROW, USERK, USERLF, SETREG, USERRG, USEREQ, ELIMEQ, 0033
C LH1405, SVDRS2, QRBD, G1,G2, DIFF, DIAREG, DIAGA, SETGA1, 0034
C SETVAL, LDPETC, LDP, NNLS, CVNEQ, FISHNI, GAMLN, 0035
C BETAIN, PLPRIN, USEROU, MOMENT, MOMOUT, RUNRES, GETPRU, GETYLY, 0036
C PGAUSS, PLRES, SETSGN, ANPEAK, UPDSGN, UPDDON, FFLAT, UPDLLS, 0037
C USERWT 0038
C----------------------------------------------------------------------- 0039
DOUBLE PRECISION PRECIS, RANGE 0040
DOUBLE PRECISION A, AA, AEQ, AINEQ, PIVOT, REG, RHSNEQ, 0041
1 S, SOLBES, SOLUTN, SSCALE, VALPCV, VALPHA, VK1Y1, WORK 0042
LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0043
1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0044
LOGICAL LBIND 0045
C 0046
C*********************************************************************** 0047
C THE INSTRUCTIONS SET OFF BY ASTERISKS DESCRIBE ALL POSSIBLE CHANGES 0048
C THAT YOU MAY HAVE TO MAKE IN THIS MAIN SUBPROGRAM. (SEE ALSO THE 0049
C CHANGES IN THE BLOCK DATA AND USER SUBPROGRAMS.) THESE CHANGES 0050
C IN THE MAIN SUBPROGRAM ARE ONLY NECESSARY IF YOU CHANGE MY, MA, 0051
C MG, MREG, MINEQ, MEQ, MDONE, OR MWORK IN THE DATA STATEMENT 0052
C BELOW. IF YOU DO, THEN THE FOLLOWING DIMENSIONS MUST BE 0053
C READJUSTED AS DESCRIBED BELOW - 0054
C 0055
COMMON /MS1/ AA
COMMON /MS2/ AINEQ
COMMON /MS3/ A
COMMON /MS4/ REG
COMMON /MS5/ AEQ
COMMON /MS6/ WORK
DIMENSION T(8192), SQRTW(8192), Y(8192), EXACT(8192), YLYFIT(8192)
DIMENSION G(503), CQUAD(503), VK1Y1(503), S(503,3), VALPHA(503),
1 VALPCV(503), SOLUTN(503), IISIGN(503), SOLBES(503),
2 AA(503,503), SSCALE(503)
DIMENSION AINEQ(501,503), RHSNEQ(501), LBIND(501)
DIMENSION A(503,503), IWORK(503)
DIMENSION REG(504,503)
DIMENSION AEQ(11,503), PIVOT(11)
DIMENSION WORK(253508)
DIMENSION LSDONE(90,3,2), VDONE(90) 0065
C 0066
C THE ABOVE DIMENSION STATEMENTS MUST BE ADJUSTED AS FOLLOWS - 0067
C 0068
C DIMENSION T(MY), SQRTW(MY), Y(MY), EXACT(MY), YLYFIT(MY) 0069
C DIMENSION G(MG), CQUAD(MG), VK1Y1(MG), S(MG,3), VALPHA(MG), 0070
C 1 VALPCV(MG), SOLUTN(MG), IISIGN(MG), SOLBES(MG), 0071
C 2 AA(MG,MG), SSCALE(MG) 0072
C DIMENSION AINEQ(MINEQ,MG), RHSNEQ(MINEQ), LBIND(MINEQ) 0073
C DIMENSION A(MA,MG), IWORK(MA) 0074
C DIMENSION REG(MREG,MG) 0075
C DIMENSION AEQ(MEQ,MG), PIVOT(MEQ) 0076
C DIMENSION WORK(MWORK) 0077
C DIMENSION LSDONE(MDONE,3,2), VDONE(MDONE) 0078
C*********************************************************************** 0079
C 0080
COMMON /DBLOCK/ PRECIS, RANGE 0081
COMMON /SBLOCK/ DFMIN, SRMIN, 0082
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 0083
2 SRANGE 0084
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0085
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0086
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0087
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0088
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0089
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0090
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0091
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0092
2 LUSER(30) 0093
C 0094
C*********************************************************************** 0095
C YOU CAN SAVE STORAGE BY MAKING THE INTEGERS IN THE FOLLOWING DATA 0096
C STATEMENT AS SMALL AS THE SIZE OF YOUR PROBLEM WILL ALLOW. (SEE 0097
C USERS MANUAL FOR THE MINIMUM ALLOWABLE VALUES.) 0098
C 0099
DATA MY/8192/, MA/503/, MG/503/, MREG/504/, MINEQ/501/, MEQ/11/,
1 MDONE/90/, MWORK/253508/
C 0102
C IF YOU CHANGE THE ABOVE DATA STATEMENT, THEN THE DIMENSION 0103
C STATEMENTS ABOVE MUST BE READJUSTED, AS DESCRIBED ABOVE. 0104
C 0105
C THIS IS THE END OF ALL POSSIBLE CHANGES THAT YOU MIGHT HAVE TO MAKE 0106
C IN THE MAIN PROGRAM, 0107
C EXCEPT THAT IF YOUR SYSTEM DOES NOT AUTOMATICALLY 0108
C OPEN INPUT AND OUTPUT FILES FOR YOU, THEN YOU MIGHT HAVE TO OPEN 0109
C THEM HERE AND GIVE THEM THE NUMBERS NIN (FOR THE INPUT) AND 0110
C NOUT (FOR THE OUTPUT), WHERE NIN AND NOUT HAVE BEEN SET IN THE 0111
C BLOCK DATA SUBPROGRAM. 0112
C IN ADDITION, IF YOU ARE GOING TO INPUT IUNIT.GE.0, THEN YOU MAY 0113
C HAVE TO OPEN A TEMPORARY SCRATCH FILE NUMBERED IUNIT. DO THIS 0114
C DIRECTLY AFTER STATEMENT 100. THIS IS NOT NECESSARY IF IUNIT 0115
C IS NEGATIVE OR IF YOUR SYSTEM OPENS FILES AUTOMATICALLY. 0116
C*********************************************************************** 0117
C 0118
C----------------------------------------------------------------------- 0119
C INITIALIZE VARIABLES 0120
C----------------------------------------------------------------------- 0121
CALL INIT 0122
C----------------------------------------------------------------------- 0123
C READ INPUT DATA 0124
C----------------------------------------------------------------------- 0125
100 CALL INPUT (EXACT,G,MA,MEQ,MG,MINEQ,MREG,MWORK,MY,SQRTW,T,Y) 0126
C----------------------------------------------------------------------- 0127
C SET UP QUADRATURE GRID 0128
C----------------------------------------------------------------------- 0129
CALL SETGRD (CQUAD,G,GMNMX,IGRID,IQUAD,MG,NG,NOUT) 0130
C----------------------------------------------------------------------- 0131
C CALCULATE SIMULATED DATA 0132
C----------------------------------------------------------------------- 0133
IF (SIMULA) CALL USERSI (EXACT,G,MG,MY,SQRTW,T,Y) 0134
C----------------------------------------------------------------------- 0135
C WRITE OUT SIMULATED DATA 0136
C----------------------------------------------------------------------- 0137
IF (SIMULA) CALL WRITYT (EXACT,G,IPRINT,IUSROU,IWT,MG,NOUT,NY, 0138
1 PRY,SIMULA,SQRTW,T,Y) 0139
C----------------------------------------------------------------------- 0140
C PUT SECOND CURVE TO BE PLOTTED WITH SOLUTION IN EXACT. 0141
C----------------------------------------------------------------------- 0142
IF (.NOT.ONLY1) CALL USERSX (EXACT,G,MG) 0143
NINEQ=0 0144
NGL=NG+NLINF 0145
NGLP1=NGL+1 0146
C----------------------------------------------------------------------- 0147
C SET SPECIAL USER-SUPPLIED INEQUALITY CONSTRAINTS 0148
C----------------------------------------------------------------------- 0149
IF (DOUSNQ) CALL USERNQ (AINEQ,MG,MINEQ) 0150
C----------------------------------------------------------------------- 0151
C SET NG NONNEGATIVITY CONSTRAINTS AT ALL NG GRID POINTS 0152
C----------------------------------------------------------------------- 0153
IF (NONNEG) CALL SETNNG (AINEQ,MINEQ,NG,NGLP1,NINEQ,NOUT) 0154
IF (IWT.EQ.1 .OR. IWT.EQ.4) GO TO 200 0155
C----------------------------------------------------------------------- 0156
C DO A COMPLETE PRELIMINARY UNWEIGHTED ANALYSIS TO GET A SMOOTH FIT 0157
C TO THE DATA. THIS SMOOTH CURVE IS THEN USED TO CALCULATE THE 0158
C WEIGHTS. 0159
C----------------------------------------------------------------------- 0160
CALL ANALYZ (1, 0161
1 A,AA,AEQ,AINEQ,CQUAD,EXACT,G,IISIGN,IWORK,LBIND,LSDONE,MA, 0162
2 MDONE,MEQ,MG,MINEQ,MREG,MWORK,MY,PIVOT,REG,RHSNEQ,S,SOLBES, 0163
3 SOLUTN,SQRTW,SSCALE,T,VALPCV,VALPHA,VDONE,VK1Y1,WORK, 0164
4 Y,YLYFIT) 0165
C----------------------------------------------------------------------- 0166
C CALCULATE SQRTW (SQUARE ROOT OF LEAST SQUARES WEIGHTS). 0167
C----------------------------------------------------------------------- 0168
CALL SETWT ( 0169
1 CQUAD,G,IUNIT,IWT,MWORK,MY,NERFIT,NG,NGL,NLINF,NOUT,NY,PRWT, 0170
2 SOLBES,SQRTW,SRANGE,SSCALE,T,WORK,Y,YLYFIT) 0171
C----------------------------------------------------------------------- 0172
C DO FINAL WEIGHTED ANALYSIS. 0173
C----------------------------------------------------------------------- 0174
200 CALL ANALYZ (2, 0175
1 A,AA,AEQ,AINEQ,CQUAD,EXACT,G,IISIGN,IWORK,LBIND,LSDONE,MA, 0176
2 MDONE,MEQ,MG,MINEQ,MREG,MWORK,MY,PIVOT,REG,RHSNEQ,S,SOLBES, 0177
3 SOLUTN,SQRTW,SSCALE,T,VALPCV,VALPHA,VDONE,VK1Y1,WORK, 0178
4 Y,YLYFIT) 0179
IF (.NOT.LAST) GO TO 100 0180
STOP 0181
END 0182
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0183
C BLOCK DATA SUBPROGRAM. 0184
BLOCK DATA 0185
DOUBLE PRECISION PRECIS, RANGE 0186
LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0187
1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0188
COMMON /DBLOCK/ PRECIS, RANGE 0189
COMMON /SBLOCK/ DFMIN, SRMIN, 0190
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 0191
2 SRANGE 0192
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0193
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0194
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0195
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0196
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0197
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0198
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0199
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0200
2 LUSER(30) 0201
C 0202
C*********************************************************************** 0203
C YOU MUST SET THE FOLLOWING 4 VARIABLES TO VALUES APPROPRIATE FOR 0204
C YOUR COMPUTER. (SEE USERS MANUAL.) 0205
C DATA RANGE/1.E35/, SRANGE/1.E35/, NIN/5/, NOUT/6/!SP 0206
DATA RANGE/1.D35/, SRANGE/1.E35/, NIN/5/, NOUT/6/ 0207
C*********************************************************************** 0208
C 0209
C 0210
C*********************************************************************** 0211
C CONTIN WILL USE THE VALUES OF THE CONTROL VARIABLES THAT ARE IN THE 0212
C DATA STATEMENTS BELOW UNLESS YOU INPUT DIFFERENT VALUES. TO SAVE 0213
C INPUTTING THESE OFTEN, YOU CAN CHANGE THE VALUES BELOW TO THOSE 0214
C THAT YOU WILL USUALLY USE. 0215
C THE FOLLOWING DATA STATEMENTS CONTAIN (IN ALPHABETICAL ORDER) THE 0216
C REAL, INTEGER, AND LOGICAL CONTROL VARIABLES, IN THAT ORDER. 0217
DATA ALPST/2*0./, DFMIN/2./, GMNMX/2*0./, PLEVEL/4*.5/, 0218
1 RSVMNX/2*1., 2*0./, RUSER/551*0./, 0219
2 SRMIN/.01/ 0220
DATA ICRIT/2*1/, 0221
1 IFORMT/1H(, 1H5, 1HE, 1H1, 1H5, 1H., 1H6, 1H), 62*1H /, 0222
2 IFORMW/1H(, 1H5, 1HE, 1H1, 1H5, 1H., 1H6, 1H), 62*1H /, 0223
3 IFORMY/1H(, 1H5, 1HE, 1H1, 1H5, 1H., 1H6, 1H), 62*1H /, 0224
4 IGRID/2/, IPLFIT/2*2/, IPLRES/2*2/, IPRINT/2*4/, IQUAD/3/, 0225
5 IUNIT/-1/, IUSER/9*0, 2, 7*0, 50, 32*0/, IUSROU/2*0/, IWT/1/, 0226
6 LINEPG/60/, LSIGN/16*0/, MIOERR/5/, MOMNMX/-1, 3/, MPKMOM/5/, 0227
7 MQPITR/35/, NENDZ/2, 2/, NEQ/0/, NERFIT/10/, NFLAT/8*0/, NG/31/, 0228
8 NINTT/1/, NLINF/0/, NNSGN/2*0/, NORDER/2/, NQPROG/6, 6/, 0229
9 NSGN/4*0/ 0230
DATA DOCHOS/.TRUE./, DOMOM/.TRUE./, DOUSIN/.TRUE./, 0231
1 DOUSNQ/.FALSE./, LAST/.TRUE./, 0232
2 LUSER/30*.FALSE./, NEWPG1/.FALSE./, NONNEG/.TRUE./, 0233
3 ONLY1/.TRUE./, PRWT/.TRUE./, PRY/.TRUE./, 0234
4 SIMULA/.FALSE./ 0235
C*********************************************************************** 0236
C 0237
C 0238
C*********************************************************************** 0239
C IF YOU MAKE ANY CHANGES TO CONTIN, YOU SHOULD PUT A NAME (UP TO 6 0240
C CHARACTERS) IN IAPACK TO UNIQUELY IDENTIFY YOUR VERSION OF 0241
C CONTIN. THIS WILL BE PRINTED IN THE HEADING OF VARIOUS PARTS OF 0242
C THE OUTPUT. YOU MUST SPECIFY THE NAME IN THE FOLLOWING STATEMENT 0243
DATA IAPACK/1H ,1HP, 1HC, 1HS, 1H-, 1H1/ 0244
C*********************************************************************** 0245
C 0246
END 0247
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0248
C SUBROUTINE USEREQ. THIS IS A USER-SUPPLIED ROUTINE (NEEDED 0249
C WHEN NEQ IS INPUT POSITIVE) TO PUT THE EQUALITY-CONSTRAINT 0250
C MATRIX IN THE FIRST NEQ ROWS AND NGL COLUMNS OF AEQ, AND 0251
C THE RIGHT-HAND-SIDE OF THE EQUALITIES IN COLUMN NGLP1=NGL+1. 0252
C CQUAD CONTAINS THE COEFFICIENTS OF THE QUADRATURE FORMULA THAT WERE 0253
C SET IN SETGRD. 0254
C NOTE - IF WEIGHTS ARE TO BE CALCULATED (I.E., IF IWT=2,3, OR 5), THEN 0255
C USEREQ WILL BE CALLED TWICE. THEREFORE IT IS BEST TO HAVE ANY 0256
C DATA NEEDED BY USEREQ READ IN ONLY ONCE AND STORED, E.G., IN 0257
C RUSER, AS ILLUSTRATED BELOW. 0258
C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 0259
C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 0260
C BELOW IS ILLUSTRATED THE CASE WHERE THE END POINTS AND THE 0261
C INTEGRAL OVER THE SOLUTION CAN BE CONSTRAINED - 0262
C IF NEQ=1, THEN THE SOLUTION IS CONSTRAINED TO BE RUSER(1) AT 0263
C THE GRID POINT NG (THE LAST POINT). 0264
C IF NEQ=2, THEN, IN ADDITION, THE SOLUTION IS CONSTRAINED TO 0265
C RUSER(2) AT GRID POINT 1 0266
C IF NEQ=3, THEN, IN ADDITION, THE INTEGRAL OVER THE SOLUTION 0267
C (USING THE QUADRATURE APPROXIMATION) IS CONSTRAINED TO BE 0268
C RUSER(6). NEQ = 1, 2 OR 3 AND RUSER(J), J=1, (AND 2 AND 6 0269
C IF NEQ=2 OR 3) WOULD HAVE TO HAVE BEEN PREVIOUSLY INPUT. 0270
SUBROUTINE USEREQ (AEQ,CQUAD,MEQ,MG) 0271
DOUBLE PRECISION PRECIS, RANGE 0272
DOUBLE PRECISION AEQ 0273
LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0274
1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0275
DIMENSION AEQ(MEQ,MG), CQUAD(MG) 0276
COMMON /DBLOCK/ PRECIS, RANGE 0277
COMMON /SBLOCK/ DFMIN, SRMIN, 0278
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 0279
2 SRANGE 0280
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0281
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0282
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0283
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0284
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0285
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0286
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0287
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0288
2 LUSER(30) 0289
C ZERO=0.E0!SP 0290
ZERO=0.D0 0291
C ONE=1.E0!SP 0292
ONE=1.D0 0293
C----------------------------------------------------------------------- 0294
C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 0295
C FOR YOUR APPLICATION. 0296
C----------------------------------------------------------------------- 0297
IF (NEQ.GE.1 .AND. NEQ.LE.3) GO TO 105 0298
5105 FORMAT (/6H NEQ =,I3,28H IS NOT 1, 2 OR 3 IN USEREQ.) 0299
WRITE (NOUT,5105) NEQ 0300
STOP 0301
105 L=MIN0(NEQ,2) 0302
DO 110 J=1,L 0303
DO 120 K=1,NGL 0304
AEQ(J,K)=ZERO 0305
120 CONTINUE 0306
AEQ(J,NGLP1)=RUSER(J) 0307
110 CONTINUE 0308
AEQ(1,NG)=ONE 0309
IF (NEQ .GT. 1) AEQ(2,1)=ONE 0310
IF (NEQ .NE. 3) GO TO 800 0311
DO 130 K=1,NGL 0312
AEQ(3,K)=ZERO 0313
IF (K .LE. NG) AEQ(3,K)=CQUAD(K) 0314
130 CONTINUE 0315
AEQ(3,NGLP1)=RUSER(6) 0316
800 RETURN 0317
END 0318
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0319
C FUNCTION USEREX. THIS IS A USER-SUPPLIED FUNCTION (ONLY USED 0320
C IF SIMULA=.TRUE.) TO EVALUATE THE NOISE-FREE VALUE OF THE 0321
C SIMULATED DATA AT T(IROW) (THE INDEPENDENT VARIABLE AT THE 0322
C DATA POINT IROW). 0323
C BELOW IS ILLUSTRATED THE CASE OF THE SOLUTION BEING COMPOSED OF 0324
C A SUM OF IUSER(11) DIRAC DELTA FUNCTIONS LOCATED AT 0325
C RUSER(25), RUSER(27), ... , RUSER(23+2*IUSER(11)) 0326
C WITH RESPECTIVE AMPLITUDES 0327
C RUSER(26), RUSER(28), ... , RUSER(24+2*IUSER(11)). 0328
C IN ADDITION, THE NLINF FUNCTIONS IN USERLF ARE ADDED WITH 0329
C RESPECTIVE AMPLITUDES 0330
C RUSER(41), RUSER(42), ... , RUSER(40+NLINF). 0331
C IUSER(11) AND NLINF CANNOT BOTH BE LESS THAN ONE. 0332
C IUSER(11) CANNOT EXCEED 7. 0333
C NLINF CANNOT EXCEED 10. 0334
C WHEN IWT=5, THE NORMALIZATION CONSTANT, RUSER(14), IS COMPUTED 0335
C FOR THE SPECIAL CASE OF PHOTON CORRELATION SPECTROSCOPY SO THAT 0336
C GAMMA = RUSER(26)+RUSER(28)+...+RUSER(24+2*IUSER(11)), WHERE 0337
C GAMMA IS GIVEN IN EQ.(6) IN MAKROMOL. CHEM., VOL. 180, P. 201 0338
C (1979). 0339
C----------------------------------------------------------------------- 0340
C CALLS SUBPROGRAMS - USERLF, USERK, ERRMES 0341
C WHICH IN TURN CALL - USERTR 0342
C----------------------------------------------------------------------- 0343
FUNCTION USEREX (IROW,T,MY,G,MG) 0344
DOUBLE PRECISION PRECIS, RANGE 0345
LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0346
1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0347
DIMENSION T(MY), IHOLER(6), ADUM(1), G(MG) 0348
COMMON /DBLOCK/ PRECIS, RANGE 0349
COMMON /SBLOCK/ DFMIN, SRMIN, 0350
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 0351
2 SRANGE 0352
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0353
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0354
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0355
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0356
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0357
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0358
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0359
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0360
2 LUSER(30) 0361
DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HE, 1HX/ 0362
C----------------------------------------------------------------------- 0363
C THE FOLLOWING STATEMENTS SHOULD BE REPLACED WITH THE ONES 0364
C APPROPRIATE FOR YOUR SIMULATION. 0365
C----------------------------------------------------------------------- 0366
IF (NLINF.GT.10 .OR. IUSER(11).GT.8 .OR. 0367
1 MAX0(NLINF,IUSER(11)).LE.0) CALL ERRMES (1,.TRUE.,IHOLER,NOUT) 0368
USEREX=0. 0369
IF (NLINF .LE. 0) GO TO 150 0370
DO 110 J=1,NLINF 0371
JJ=J 0372
IF (ABS(RUSER(J+40)) .GT. 0.) USEREX=USEREX+RUSER(J+40)* 0373
1 USERLF(IROW,JJ,T,MY) 0374
110 CONTINUE 0375
150 IF (IUSER(11) .LE. 0) GO TO 800 0376
JJ=IUSER(11) 0377
IF (IROW .GT. 1) GO TO 158 0378
RUSER(14)=1. 0379
IF (IWT .NE. 5) GO TO 158 0380
RUSER(14)=0. 0381
AMPSUM=0. 0382
DO 155 J=1,JJ 0383
AMPSUM=AMPSUM+RUSER(2*J+24) 0384
ADUM(1)=0. 0385
G(NG+1)=RUSER(2*J+23) 0386
RUSER(14)=RUSER(14)+RUSER(2*J+24)*USERK(1,ADUM,NG+1,G) 0387
155 CONTINUE 0388
IF (RUSER(14) .LE. 0.) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) 0389
RUSER(14)=AMPSUM/RUSER(14) 0390
158 DO 160 J=1,JJ 0391
G(NG+1)=RUSER(2*J+23) 0392
USEREX=USEREX+RUSER(2*J+24)*USERK(IROW,T,NG+1,G)*RUSER(14) 0393
160 CONTINUE 0394
800 RETURN 0395
END 0396
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0397
C SUBROUTINE USERGR. THIS IS A USER-SUPPLIED ROUTINE (ONLY CALLED 0398
C WHEN IGRID IS INPUT AS 3) FOR COMPUTING A NONSTANDARD GRID G 0399
C AND THE QUADRATURE WEIGHTS CQUAD. 0400
C IN THE EXAMPLE BELOW, G IS SIMPLY READ IN. CQUAD IS 0401
C THEN COMPUTED FOR THE TRAPEZOIDAL RULE. 0402
C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 0403
C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 0404
C----------------------------------------------------------------------- 0405
C CALLS SUBPROGRAMS - CQTRAP 0406
C----------------------------------------------------------------------- 0407
SUBROUTINE USERGR (G,CQUAD,MG) 0408
DOUBLE PRECISION PRECIS, RANGE 0409
LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0410
1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0411
DIMENSION G(MG), CQUAD(MG) 0412
COMMON /DBLOCK/ PRECIS, RANGE 0413
COMMON /SBLOCK/ DFMIN, SRMIN, 0414
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 0415
2 SRANGE 0416
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0417
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0418
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0419
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0420
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0421
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0422
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0423
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0424
2 LUSER(30) 0425
C----------------------------------------------------------------------- 0426
C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 0427
C FOR YOUR APPLICATION. 0428
C----------------------------------------------------------------------- 0429
READ (NIN,5100) (G(J),J=1,NG) 0430
5100 FORMAT (5E15.6) 0431
C----------------------------------------------------------------------- 0432
C CQTRAP USES G TO PUT TRAPEZOIDAL RULE WEIGHTS IN CQUAD. 0433
C----------------------------------------------------------------------- 0434
CALL CQTRAP (G,CQUAD,NG) 0435
RETURN 0436
END 0437
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0438
C SUBROUTINE USERIN. THIS IS A USER-SUPPLIED ROUTINE (ONLY CALLED 0439
C WHEN DOUSIN=.TRUE.) THAT IS CALLED RIGHT AFTER THE INITIALIZATION 0440
C AND INPUT OF THE COMMON VARIABLES, T, Y, AND THE LEAST-SQUARES 0441
C WEIGHTS. THEREFORE, IT CAN BE USED TO MODIFY ANY OF THESE. 0442
C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 0443
C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 0444
C BELOW IS ILLUSTRATED TWO PREPARATION OF INPUT DATA FOR KERNELS OF 0445
C THE GENERAL FORM - 0446
C 0447
C USERK(T,G)=FORMF2(G)*G**RUSER(23)*EXP(-RUSER(21)*T*G**RUSER(22)) 0448
C 0449
C WHERE RUSER(21) AND RUSER(22) ARE NOT ZERO. THIS INCLUDES 0450
C LAPLACE TRANSFORMS AND MANY APPLICATIONS IN PHOTON CORRELATION 0451
C SPECTROSCOPY AS SPECIAL CASES. 0452
C 0453
C THE FOLLOWING ARE USED INTERNALLY BY USERIN AND USERK, 0454
C AND CANNOT BE USED BY YOU FOR ANY OTHER PURPOSES - 0455
C IUSER(10), IUSER(18), LUSER(3), LUSER(4), 0456
C RUSER(J), J = 10,15,...,50+NG (AND THEREFORE NG 0457
C CANNOT EXCEED 50). 0458
C 0459
C THE FOLLOWING ARE THE NECESSARY INPUT - 0460
C 0461
C DOUSIN = T (DOUSIN MUST ALWAYS BE .TRUE..) 0462
C 0463
C LUSER(3) = T, TO HAVE FORMF2, THE SQUARED FORM FACTORS, COMPUTED IN 0464
C USERK. 0465
C = F, TO SET ALL THE FORMF2 TO 1. (AS WOULD BE APPROPRIATE 0466
C WITH LAPLACE TRANSFORMS). 0467
C RUSER(24) MAY BE NECESSARY INPUT TO SPECIFY THE FORM FACTOR (E.G., 0468
C THE WALL THICKNESS OF A HOLLOW SPHERE) IF LUSER(3)=T. SEE 0469
C COMMENTS IN USERK. 0470
C IUSER(18) MAY BE NECESSARY INPUT IF LUSER(3)=T (E.G., TO SPECIFY THE 0471
C NUMBER OF POINTS OVER WHICH THE SQUARED FORM FACTORS WILL 0472
C BE AVERAGED). SEE COMMENTS IN USERK. 0473
C 0474
C RUSER(16) = WAVELENGTH OF INCIDENT LIGHT (IN NANOMETERS), 0475
C = 0, IF RUSER(20), THE MAGNITUDE OF THE SCATTERING VECTOR 0476
C (IN CM**-1), IS NOT TO BE COMPUTED. WHEN 0477
C RUSER(16)=0, RUSER(15) AND RUSER(17) NEED NOT BE 0478
C INPUT, AND CONTIN WILL SET RUSER(21)=1 0479
C (AS APPROPRIATE WITH LAPLACE TRANSFORMS). 0480
C 0481
C RUSER(15) = REFRACTIVE INDEX. 0482
C RUSER(17) = SCATTERING ANGLE (IN DEGREES). 0483
C 0484
C 0485
C IUSER(10) SELECTS SPECIAL CASES OF USERK FOR MORE CONVENIENT USE. 0486
C 0487
C IUSER(10) = 1, FOR MOLECULAR WEIGHT DISTRIBUTIONS FROM PCS 0488
C (WHERE THE SOLUTION, S(G), IS SUCH THAT S(G)DG IS 0489
C THE WEIGHT FRACTION WITH MOLECULAR WEIGHT BETWEEN 0490
C G AND G+DG). 0491
C CONTIN SETS - 0492
C RUSER(23) = 1., 0493
C RUSER(21) = RUSER(18)*RUSER(20)**2. 0494
C (SEE ABOVE DISCUSSION OF RUSER(16).) 0495
C YOU MUST INPUT - 0496
C RUSER(18) TO SATISFY THE EQUATION (IN CGS UNITS) - 0497
C (DIFFUSION COEFF.)=RUSER(18)*(MOL. WT.)**RUSER(22). 0498
C RUSER(22) (MUST ALSO BE INPUT, TYPICALLY ABOUT -.5). 0499
C 0500
C IUSER(10) = 2, FOR DIFFUSION-COEFFICIENT DISTRIBUTONS OR LAPLACE 0501
C TRANSFORMS (WHERE G IS DIFF. COEFF. IN CM**2/SEC 0502
C OR, E.G., TIME CONSTANT). 0503
C CONTIN SETS - 0504
C RUSER(23) = 0., 0505
C RUSER(22) = 1., 0506
C RUSER(21) = RUSER(20)**2 (SEE ABOVE DISCUSSION 0507
C OF RUSER(16).). 0508
C 0509
C IUSER(10) = 3, FOR SPHERICAL-RADIUS DISTRIBUTIONS, ASSUMING THE 0510
C EINSTEIN-STOKES RELATION (WHERE THE SOLUTION, S(G), 0511
C IS SUCH THAT S(G)DG IS THE WEIGHT FRACTION OF 0512
C PARTICLES WITH RADIUS (IN CM) BETWEEN G AND G+DG. 0513
C WEIGHT-FRACTION DISTRIBUTIONS YIELD BETTER SCALED 0514
C PROBLEMS THAN NUMBER-FRACTION DISTRIBUTIONS, WHICH 0515
C WOULD REQUIRE RUSER(23)=6.) 0516
C CONTIN SETS - 0517
C RUSER(23) = 3., 0518
C RUSER(22) = -1., 0519
C RUSER(21) = RUSER(20)**2*(BOLTZMANN CONST.)* 0520
C RUSER(18)/(.06*PI*RUSER(19)). 0521
C (SEE ABOVE DISCUSSION OF RUSER(16).) 0522
C YOU MUST HAVE INPUT - 0523
C RUSER(18) = TEMPERATURE (IN DEGREES KELVIN), 0524
C RUSER(19) = VISCOSITY (IN CENTIPOISE). 0525
C 0526
C IUSER(10) = 4, FOR GENERAL CASE, WHERE YOU MUST HAVE INPUT - 0527
C RUSER(J), J = 21, 22, 23. 0528
C 0529
C 0530
C RUSER(J) FOR J = 18, 19, 21, 22, 23 - SEE THE ABOVE INSTRUCTIONS 0531
C FOR THE VALUE OF IUSER(10) 0532
C THAT YOU ARE USING. 0533
C 0534
C 0535
C RUSER(10) SPECIFIES HOW THE INPUT Y VALUES WILL BE CONVERTED (E.G., 0536
C TO THE NORMALIZED FIRST-ORDER CORRELATION FUNCTION). 0537
C RUSER(10) = 0., FOR NO CONVERSION (I.E., IF THE INPUT Y VALUES 0538
C ARE ALREADY IN THE REQUIRED FORM, AS WITH 0539
C LAPLACE TRANSFORMS). 0540
C RUSER(10) NEGATIVE, WHEN THE INPUT Y ARE NORMALIZED 2ND-ORDER 0541
C CORRELATION FUNCTIONS MINUS 1 (AS WITH TEST 0542
C DATA SET 1) THIS CAUSES - 0543
C Y=SIGN(SQRT(ABS(Y)),Y). 0544
C (SEE THE USERS MANUAL FOR AN 0545
C EXPLANATION OF WHY SIGN() IS USED.) 0546
C RUSER(10) POSITIVE, WHEN THE INPUT Y ARE 2ND-ORDER CORRELATION 0547
C FUNCTIONS. RUSER(10) = THE BACKGROUND TERM 0548
C (E.G., RUSER(10)=1 WHEN THE INPUT Y ARE 0549
C NORMALIZED 2ND-ORDER CORREL. FCTNS.). THIS 0550
C CAUSES - 0551
C Y=SIGN(SQRT(ABS(X)),X), WHERE 0552
C X=Y/RUSER(10)-1. 0553
C 0554
C----------------------------------------------------------------------- 0555
C CALLS SUBPROGRAMS - ERRMES 0556
C----------------------------------------------------------------------- 0557
SUBROUTINE USERIN (T,Y,SQRTW,MY) 0558
DOUBLE PRECISION PRECIS, RANGE 0559
LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0560
1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0561
DIMENSION T(MY), Y(MY), SQRTW(MY) 0562
DIMENSION IHOLER(6) 0563
COMMON /DBLOCK/ PRECIS, RANGE 0564
COMMON /SBLOCK/ DFMIN, SRMIN, 0565
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 0566
2 SRANGE 0567
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0568
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0569
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0570
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0571
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0572
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0573
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0574
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0575
2 LUSER(30) 0576
DATA IHOLER /1HU, 1HS, 1HE, 1HR, 1HI, 1HN/ 0577
C----------------------------------------------------------------------- 0578
C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 0579
C FOR YOUR APPLICATION. 0580
C----------------------------------------------------------------------- 0581
C INITIALIZE FLAG FOR CALCULATION OF FORM FACTORS IN USERK. 0582
C----------------------------------------------------------------------- 0583
LUSER(4)=.FALSE. 0584
IF (IUSER(10).LT.1 .OR. IUSER(10).GT.4) CALL ERRMES (1,.TRUE., 0585
1 IHOLER,NOUT) 0586
IF (ABS(RUSER(10)).LE.0. .OR. SIMULA) GO TO 120 0587
C----------------------------------------------------------------------- 0588
C TRANSFORM INPUT Y VALUES TO NORMALIZED FIRST-ORDER CORRELATION FCTNS. 0589
C----------------------------------------------------------------------- 0590
DO 110 J=1,NY 0591
DUM=Y(J) 0592
IF (RUSER(10) .GT. 0.) DUM=DUM/RUSER(10)-1. 0593
Y(J)=SIGN(SQRT(ABS(DUM)),DUM) 0594
110 CONTINUE 0595
120 IF (IUSER(10) .EQ. 2) RUSER(22)=1. 0596
IF (IUSER(10) .EQ. 3) RUSER(22)=-1. 0597
C----------------------------------------------------------------------- 0598
C COMPUTE THE CONSTANTS RUSER(J), J=20,21 FOR USE IN USERK. 0599
C----------------------------------------------------------------------- 0600
IF (IUSER(10) .NE. 4) RUSER(21)=1. 0601
IF (RUSER(16) .LE. 0.) GO TO 800 0602
C----------------------------------------------------------------------- 0603
C 1.256...E8 = 4.E7*PI 8.726...E-3 = .5 RADIAN/DEGREE 0604
C----------------------------------------------------------------------- 0605
RUSER(20)=1.256637E8*RUSER(15)*SIN(8.726646E-3*RUSER(17))/ 0606
1 RUSER(16) 0607
IF (IUSER(10) .EQ. 4) GO TO 800 0608
RUSER(21)=RUSER(20)**2 0609
IF (IUSER(10) .EQ. 1) RUSER(21)=RUSER(21)*RUSER(18) 0610
IF (IUSER(10) .NE. 3) GO TO 800 0611
IF (RUSER(19) .LE. 0.) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) 0612
C----------------------------------------------------------------------- 0613
C 7.323...E-16 = (BOLTZMANN CONSTANT)/(.06*PI) 0614
C----------------------------------------------------------------------- 0615
RUSER(21)=RUSER(21)*7.323642E-16*RUSER(18)/RUSER(19) 0616
800 RETURN 0617
END 0618
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0619
C FUNCTION USERK. THIS IS A USER-SUPPLIED ROUTINE (ALWAYS NEEDED) 0620
C TO COMPUTE THE FREDHOLM KERNEL, USERK, WHICH DEPENDS ON T(JT) 0621
C (THE INDEPENDENT VARIABLE AT DATA POINT JT) AND G(JG) (THE 0622
C VALUE OF THE JG TH GRID POINT. 0623
C BELOW IS ILLUSTRATED THE CASE OF A GENERAL TYPE OF KERNEL 0624
C APPLICABLE TO PHOTON CORRELATION SPECTROSCOPY AND LAPLACE 0625
C TRANSFORMS - 0626
C 0627
C USERK(T,G)=FORMF2(G)*G**RUSER(23)*EXP(-RUSER(21)*T*G**RUSER(22)) 0628
C 0629
C SEE COMMENTS IN USERIN. 0630
C THE MEAN SQUARED FORM FACTOR, FORMF2(G), COMPUTED BELOW IS FOR THE 0631
C RAYLEIGH-DEBYE APPROXIMATION FOR HOLLOW SPHERES. 0632
C SEE COMMENTS BELOW 0633
C ON HOW YOU CAN MODIFY THIS FOR ANOTHER FORM FACTOR. 0634
C RUSER(24) = THICKNESS OF THE WALLS OF THE HOLLOW SPHERES (IN CM). 0635
C = 0 FOR FULL SPHERE (I.E., FOR WALL THICKNESS = RADIUS OF 0636
C SPHERE. 0637
C IUSER(18) DETERMINES THE NUMBER OF POINTS OVER WHICH THE FORM 0638
C FACTOR WILL BE AVERAGED, AS EXPLAINED IN THE COMMENTS BELOW. 0639
C----------------------------------------------------------------------- 0640
C CALLS SUBPROGRAMS - ERRMES, USERTR 0641
C----------------------------------------------------------------------- 0642
FUNCTION USERK (JT,T,JG,G) 0643
DOUBLE PRECISION PRECIS, RANGE 0644
LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0645
1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0646
DIMENSION T(JT), G(*) 0647
DIMENSION IHOLER(6) 0648
COMMON /DBLOCK/ PRECIS, RANGE 0649
COMMON /SBLOCK/ DFMIN, SRMIN, 0650
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 0651
2 SRANGE 0652
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0653
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0654
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0655
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0656
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0657
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0658
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0659
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0660
2 LUSER(30) 0661
DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HK, 1H / 0662
IF (JT.GT.NY .OR. JG.GT.NG+1 .OR. MIN0(JT,JG).LE.0) CALL 0663
1 ERRMES (1,.TRUE.,IHOLER,NOUT) 0664
C----------------------------------------------------------------------- 0665
C THE FOLLOWING STATEMENTS SHOULD BE REPLACED BY THOSE 0666
C APPROPRIATE FOR YOUR KERNEL. 0667
C FOR EXAMPLE, FOR THE LAPLACE INTEGRAL EQUATION, YOU WOULD 0668
C SIMPLY REPLACE ALL BUT THE LAST TWO STATEMENTS BELOW BY - 0669
C USERK=EXP(-T(JT)*G(JG)) 0670
C IT MAY NOT BE NECESSARY TO GUARD AGAINST UNDERFLOW IN EXP AS IS 0671
C DONE BELOW, BUT A FEW COMPILERS ABORT AT UNDERFLOW IN EXP. 0672
C (EXMAX IS SET IN INIT TO ALOG(SRANGE), I.E., IT IS THE 0673
C LARGEST REASONABLE EXPONENT IN EXP.) 0674
C----------------------------------------------------------------------- 0675
IF (LUSER(4)) GO TO 150 0676
LUSER(4)=.TRUE. 0677
IF (.NOT.LUSER(3)) GO TO 150
C----------------------------------------------------------------------- 0678
C PUT MEAN SQUARED FORM FACTORS IN RUSER(J), J=51,...,50+NG. 0679
C IF THE FORM FACTORS OSCILLATE SO RAPIDLY THAT THERE ARE OSCILLATIONS 0680
C WITH PERIOD LESS THAN 3 OR 4 GRID POINTS, THEN YOU SHOULD USE THE 0681
C MEAN SQUARED VALUE OVER AN INTERVAL CENTERED AT EACH GRID POINT 0682
C AND EXTENDING HALFWAY TOWARD EACH OF ITS NEIGHBORS. 0683
C IUSER(18) = THE NUMBER OF POINTS ON EACH SIDE OF THE GRID POINT 0684
C OVER WHICH THE AVERAGE WILL BE TAKEN. I.E., THERE WILL BE A 0685
C TOTAL OF 2*IUSER(18)+1 POINTS IN THE AVERAGE. 0686
C IUSER(18) = 0 FOR NO AVERAGING. I.E., THE FORM FACTOR AT THE GRID 0687
C POINT WILL BE USED. 0688
C IUSER(18) = 50 IS USUALLY ADEQUATE. 0689
C NOTE THAT, IF IGRID=2 AS USUAL, THEN THE POINTS (AS WITH THE GRID) 0690
C WILL BE TAKEN IN EQUAL INTERVALS OF H(G) IN USERTR (E.G., 0691
C USUALLY IN EQUAL INTERVALS OF LOG(G)). 0692
C LUSER(3) = T WILL USE THE RAYLEIGH-DEBYE APPROXIMATION FOR 0693
C HOLLOW SPHERES FILLED WITH SOLVENT, 0694
C = F WILL SET ALL FORM FACTORS TO 1. THIS ALSO HAPPENS IF 0695
C RUSER(16).LE.0 (SINCE NO SCATTERING VECTOR WAS PUT IN 0696
C RUSER(20)). 0697
C RUSER(24) = THICKNESS OF THE WALLS OF THE HOLLOW SPHERES (IN CM). 0698
C = 0 FOR FULL SPHERE (I.E., FOR WALL THICKNESS = RADIUS OF 0699
C SPHERE. 0700
C YOU CAN CHANGE THE COMPUTATION OF RUSER BELOW TO COMPUTE THE MEAN 0701
C SQUARED FORM FACTORS THAT ARE APPROPRIATE FOR YOUR APPLICATION. 0702
C YOU CAN ALSO READ IN THE MEAN SQUARED FORM FACTORS WITH - 0703
C READ (NIN,...) (RUSER(J+50),J=1,NG) 0704
C THIS INPUT DATA WOULD HAVE TO BE BETWEEN CARD SETS 13 AND 14A 0705
C (SEE THE USERS MANUAL.) 0706
C----------------------------------------------------------------------- 0707
IF (NG .GT. 500) CALL ERRMES (2,.TRUE.,IHOLER,NOUT)
DO 120 J=1,NG 0709
IF (LUSER(3) .AND. RUSER(16).GT.0.) GO TO 121 0710
RUSER(J+50)=1. 0711
GO TO 120 0712
C----------------------------------------------------------------------- 0713
C COMPUTE AVERAGE FORM FACTOR. 0714
C----------------------------------------------------------------------- 0715
121 DELTA=0. 0716
TRSTRT=USERTR(G(J),1) 0717
TREND=TRSTRT 0718
NPTS=IUSER(18)+1 0719
IF (J.NE.1 .AND. J.NE.NG) NPTS=NPTS+IUSER(18) 0720
IF (NPTS .LE. 1) GO TO 122 0721
IF (J .NE. 1) TRSTRT=.5*(USERTR(G(J-1),1)+USERTR(G(J),1)) 0722
IF (J .NE. NG) TREND=.5*(USERTR(G(J+1),1)+USERTR(G(J),1)) 0723
DELTA=(TREND-TRSTRT)/FLOAT(NPTS-1) 0724
122 SUM=0. 0725
TR=TRSTRT-DELTA 0726
DO 125 K=1,NPTS 0727
TR=TR+DELTA 0728
GPOINT=USERTR(TR,2) 0729
C----------------------------------------------------------------------- 0730
C YOU MUST REPLACE THE STATEMENTS BETWEEN HERE AND NO. 135 WITH THOSE 0731
C APPROPRIATE FOR YOUR FORM FACTOR IF IT IS NOT THE RAYLEIGH-DEBYE 0732
C APPROXIMATION FOR A HOLLOW SPHERE. 0733
C PINNER,POUTER,PAVG = (MAGNITUDE OF SCATTERING VECTOR)*(INNER,OUTER, 0734
C AVERAGE RADIUS) 0735
C RUSER(20) = MAGNITUDE OF SCATTERING VECTOR IN CM**(-1), 0736
C GPOINT = OUTER RADIUS IN CM. 0737
C TERM = FORM FACTOR FOR G=GPOINT. 0738
C----------------------------------------------------------------------- 0739
PINNER=0. 0740
IF (RUSER(24).GT.0. .AND. RUSER(24).LT.GPOINT) PINNER= 0741
1 RUSER(20)*(GPOINT-RUSER(24)) 0742
POUTER=RUSER(20)*GPOINT 0743
PAVG=.5*(POUTER+PINNER) 0744
PDELTA=PAVG-PINNER 0745
PD2=PDELTA**2 0746
COSPAV=COS(PAVG) 0747
IF (PDELTA .LE. .2) TERM=COSPAV*PD2*(1.-PD2*(28.-PD2)/280.) 0748
1 +PAVG*SIN(PAVG)*(3.-PD2*(.5-.025*PD2)) 0749
IF (PDELTA .GT. .2) TERM=3.*(SIN(PDELTA)*(PAVG*SIN(PAVG)+ 0750
1 COSPAV)/PDELTA-COSPAV*COS(PDELTA)) 0751
TERM=TERM/(POUTER*(POUTER+PINNER)+PINNER**2) 0752
135 SUM=SUM+TERM**2 0753
125 CONTINUE 0754
RUSER(J+50)=SUM/FLOAT(NPTS) 0755
120 CONTINUE 0756
K=NG+50 0757
5120 FORMAT (/21H SQUARED FORM FACTORS/(1P10E13.3)) 0758
IF (LUSER(3) .AND. RUSER(16).GT.0.) WRITE (NOUT,5120) 0759
1 (RUSER(J),J=51,K) 0760
C----------------------------------------------------------------------- 0761
C END OF CALCULATION OF MEAN SQUARED FORM FACTORS. 0762
C----------------------------------------------------------------------- 0763
C PUT MEAN SQUARED FORM FACTOR IN FORMF2. NORMALLY THIS IS JUST 0764
C RUSER(JG+50). HOWEVER, WITH CALLS FROM USEREX, G(NG+1) CAN HAVE 0765
C ANY VALUE. IN THIS CASE, THE FORM FACTOR IS ESTIMATED BY LINEAR 0766
C INTERPOLATION USING THE TWO GRID POINTS NEAREST G(NG+1). 0767
C ADVANTAGE IS TAKEN OF THE FACT THAT G(J) IS STRICTLY MONOTONIC 0768
C FOR J=1,...,NG. 0769
C----------------------------------------------------------------------- 0770
150 FORMF2=1. 0771
IF (.NOT.LUSER(3) .OR. RUSER(16).LE.0.) GO TO 158 0772
IF (JG .EQ. NG+1) GO TO 152 0773
FORMF2=RUSER(JG+50) 0774
GO TO 158 0775
152 IF ((G(NG+1)-G(1))*(G(NG+1)-G(NG)) .GT. 0.) CALL ERRMES (3, 0776
1 .TRUE.,IHOLER,NOUT) 0777
SGN=SIGN(1.,G(2)-G(1)) 0778
DO 155 J=2,NG 0779
JJ=J 0780
IF ((G(JJ)-G(NG+1))*SGN .GE. 0.) GO TO 157 0781
155 CONTINUE 0782
157 FORMF2=RUSER(JJ+49)+(RUSER(JJ+50)-RUSER(JJ+49))*(G(NG+1)-G(JJ-1))/ 0783
1 (G(JJ)-G(JJ-1)) 0784
158 IF (AMIN1(ABS(RUSER(21)),ABS(RUSER(22))) .LE. 0.) CALL ERRMES (4, 0785
1 .TRUE.,IHOLER,NOUT) 0786
USERK=0. 0787
IF (IUSER(10) .NE. 1) GO TO 200 0788
C----------------------------------------------------------------------- 0789
C FOR MOLECULAR WEIGHT DISTRIBUTIONS. 0790
C----------------------------------------------------------------------- 0791
IF ((RUSER(22).LE.0. .AND. G(JG).LE.0.) .OR. G(JG).LT.0.) 0792
1 CALL ERRMES (5,.TRUE.,IHOLER,NOUT) 0793
IF (ABS(.5+RUSER(22)) .GT. 1.E-5) GO TO 160 0794
EX=RUSER(21)*T(JT)/SQRT(G(JG)) 0795
GO TO 170 0796
160 EX=RUSER(21)*T(JT)*G(JG)**RUSER(22) 0797
170 PREEXP=FORMF2*G(JG) 0798
GO TO 500 0799
200 IF (IUSER(10) .NE. 2) GO TO 300 0800
C----------------------------------------------------------------------- 0801
C FOR DIFFUSION-COEFFICIENT DISTRIBUTIONS OR LAPLACE TRANSFORMS. 0802
C----------------------------------------------------------------------- 0803
EX=RUSER(21)*G(JG)*T(JT) 0804
PREEXP=FORMF2 0805
GO TO 500 0806
300 IF (IUSER(10) .NE. 3) GO TO 400 0807
C----------------------------------------------------------------------- 0808
C FOR SPHERICAL-RADIUS DISTRIBUTIONS. 0809
C----------------------------------------------------------------------- 0810
IF (G(JG) .LE. 0.) CALL ERRMES (6,.TRUE.,IHOLER,NOUT) 0811
EX=RUSER(21)*T(JT)/G(JG) 0812
PREEXP=FORMF2*G(JG)**3 0813
GO TO 500 0814
C----------------------------------------------------------------------- 0815
C GENERAL FORM OF KERNEL - WHEN IUSER(10)=4. 0816
C----------------------------------------------------------------------- 0817
400 EX=0. 0818
IF (G(JG) .LE. 0.) GO TO 410 0819
C----------------------------------------------------------------------- 0820
C G(JG) IS POSITIVE. 0821
C----------------------------------------------------------------------- 0822
EX=RUSER(21)*T(JT)*G(JG)**RUSER(22) 0823
PREEXP=FORMF2*G(JG)**RUSER(23) 0824
GO TO 500 0825
410 IF (ABS(G(JG)) .LE. 0.) GO TO 420 0826
C----------------------------------------------------------------------- 0827
C G(JG) IS NEGATIVE. 0828
C----------------------------------------------------------------------- 0829
J=INT(RUSER(22)+SIGN(.5,RUSER(22))) 0830
JJ=INT(RUSER(23)+SIGN(.5,RUSER(23))) 0831
IF (AMAX1(ABS(RUSER(22)-FLOAT(J)),ABS(RUSER(23)-FLOAT(JJ))) .GT. 0832
1 1.E-5) CALL ERRMES (7,.TRUE.,IHOLER,NOUT) 0833
EX=RUSER(21)*T(JT)*G(JG)**J 0834
PREEXP=FORMF2*G(JG)**JJ 0835
GO TO 500 0836
C----------------------------------------------------------------------- 0837
C G(JG)=0. 0838
C----------------------------------------------------------------------- 0839
420 IF (RUSER(22).LT.0. .OR. RUSER(23).GT.0.) RETURN 0840
IF (RUSER(23) .LT. 0.) CALL ERRMES (8,.TRUE.,IHOLER,NOUT) 0841
PREEXP=FORMF2 0842
500 IF (EX .LT. EXMAX) USERK=PREEXP*EXP(-EX) 0843
RETURN 0844
END 0845
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0846
C FUNCTION USERLF. THIS IS A USER-SUPPLIED ROUTINE (NEEDED IF 0847
C NLINF IS INPUT POSITIVE) TO CALCULATE THE NLINF ADDITIONAL 0848
C KNOWN FUNCTIONS THAT ARE PRESENT IN UNKNOWN AMOUNTS 0849
C IN THE DATA. 0850
C JY = THE SUBSCRIPT OF THE DATA POINT. 0851
C JLINF = INDEX TELLING WHICH OF THE NLINF FUNCTIONS IS TO BE 0852
C EVALUATED. (1 .LE. JLINF .LE. NLINF) 0853
C NYDIM = THE DIMENSION OF THE ARRAY T. IT CAN BE MY OR NY. 0854
C BELOW IS ILLUSTRATED THE CASE WHERE 2 SETS OF DATA CAN BE 0855
C COMBINED AND EACH CAN HAVE A DIFFERENT CONSTANT ADDITIVE 0856
C BACKGROUND (BASELINE). THE FIRST IUSER(2) LOCATIONS IN 0857
C Y, T, ETC. ARE OCCUPIED BY THE FIRST DATA SET. THE SECOND SET 0858
C OCCUPIES LOCATIONS IUSER(2)+1 THRU NY. THEREFORE JLINF=1 OR 2, 0859
C AND USERLF=1. OR 0. DEPENDING ON WHETHER THE (JY)TH DATA POINT 0860
C BELONGS TO DATA SET JLINF OR NOT. 0861
C IF, HOWEVER, IUSER(2)=0 AND NLINF=1, THEN IT IS ASSUMED THAT 0862
C THERE IS ONLY ONE SET OF DATA AND ONLY ONE BASELINE. (THIS 0863
C SAVES YOU THE TROUBLE OF SETTING IUSER(2)=NY EVERY TIME NY 0864
C CHANGES.) 0865
C IF NLINF=0, THEN THERE WILL BE NO BASELINE. 0866
C----------------------------------------------------------------------- 0867
C CALLS SUBPROGRAMS - ERRMES 0868
C----------------------------------------------------------------------- 0869
FUNCTION USERLF (JY,JLINF,T,NYDIM) 0870
DOUBLE PRECISION PRECIS, RANGE 0871
LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0872
1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0873
DIMENSION T(NYDIM), IHOLER(6) 0874
COMMON /DBLOCK/ PRECIS, RANGE 0875
COMMON /SBLOCK/ DFMIN, SRMIN, 0876
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 0877
2 SRANGE 0878
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0879
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0880
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0881
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0882
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0883
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0884
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0885
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0886
2 LUSER(30) 0887
DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HL, 1HF/ 0888
IF (JY.GT.NY .OR. JY.LE.0) CALL ERRMES (1,.TRUE.,IHOLER,NOUT) 0889
C----------------------------------------------------------------------- 0890
C THE FOLLOWING STATEMENTS MUST BE REPLACED BY THE APPROPRIATE 0891
C ONES FOR YOUR ADDITIVE FUNCTIONS. 0892
C----------------------------------------------------------------------- 0893
IF (JLINF.LT.1 .OR. JLINF.GT.2) CALL ERRMES (2,.TRUE.,IHOLER, 0894
1 NOUT) 0895
USERLF=0. 0896
IF ((JY.LE.IUSER(2) .AND. JLINF.EQ.1) .OR. (JY.GT.IUSER(2) .AND. 0897
1 JLINF.EQ.2) .OR. IUSER(2).LE.0) USERLF=1. 0898
RETURN 0899
END 0900
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0901
C SUBROUTINE USERNQ. THIS IS A USER-SUPPLIED ROUTINE (ONLY 0902
C CALLED WHEN DOUSNQ IS INPUT AS .TRUE.) TO SET NINEQ 0903
C (THE NO. OF INEQUALITY CONSTRAINTS) AND TO PUT THE 0904
C INEQUALITY-CONSTRAINT MATRIX IN THE FIRST NINEQ 0905
C ROWS OF AINEQ AND TO PUT THE RIGHT HAND SIDES OF THE 0906
C INEQUALITIES IN COLUMN NGLP1 OF AINEQ. 0907
C I.E., (SUM FROM J=1 TO NGL OF (AINEQ(I,J)*SOLUTION(J))) .GE. 0908
C AINEQ(I,NGLP1), I=1,NINEQ. SEE EQ. (3.6). 0909
C NOTE - IF THE SOLUTION IS TO BE NONNEGATIVE AT ALL NG GRID 0910
C POINTS, DO NOT USE USERNQ TO SET THESE CONSTRAINTS - 0911
C SIMPLY INPUT NONNEG=.TRUE. INSTEAD. 0912
C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 0913
C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 0914
C BELOW IS ILLUSTRATED THE CASE WHERE THE NLINF COEFFICIENTS OF 0915
C THE NLINF LINEAR FUNCTIONS ARE CONSTRAINED TO BE 0916
C NONNEGATIVE. 0917
C----------------------------------------------------------------------- 0918
C CALLS SUBPROGRAMS - ERRMES 0919
C----------------------------------------------------------------------- 0920
SUBROUTINE USERNQ (AINEQ,MG,MINEQ) 0921
DOUBLE PRECISION PRECIS, RANGE 0922
DOUBLE PRECISION AINEQ, ONE, ZERO 0923
LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 0924
1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0925
DIMENSION AINEQ(MINEQ,MG) 0926
DIMENSION IHOLER(6) 0927
COMMON /DBLOCK/ PRECIS, RANGE 0928
COMMON /SBLOCK/ DFMIN, SRMIN, 0929
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 0930
2 SRANGE 0931
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 0932
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 0933
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 0934
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 0935
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 0936
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 0937
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 0938
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 0939
2 LUSER(30) 0940
DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HN, 1HQ/ 0941
C ZERO=0.E0!SP 0942
ZERO=0.D0 0943
C ONE=1.E0!SP 0944
ONE=1.D0 0945
C----------------------------------------------------------------------- 0946
C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 0947
C FOR YOUR APPLICATION. 0948
C----------------------------------------------------------------------- 0949
IF (NLINF .LE. 0) RETURN 0950
NINEQ=NLINF 0951
IF (NINEQ .GT. MINEQ) CALL ERRMES (1,.TRUE.,IHOLER,NOUT) 0952
DO 110 J=1,NINEQ 0953
DO 120 K=1,NGLP1 0954
AINEQ(J,K)=ZERO 0955
120 CONTINUE 0956
K=NG+J 0957
AINEQ(J,K)=ONE 0958
110 CONTINUE 0959
RETURN 0960
END 0961
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0962
C SUBROUTINE USEROU. THIS IS A USER-SUPPLIED ROUTINE (ONLY CALLED 0963
C WHEN DOUSOU=.TRUE.) THAT YOU CAN USE TO PRODUCE YOUR OWN EXTRA 0964
C OUTPUT. 0965
C G CONTAINS THE NG GRID POINTS. 0966
C SOL CONTAINS NG+NLINF VALUES - FIRST THE NG SOLUTION VALUES AND THEN 0967
C THE NLINF LINEAR COEFFICIENTS. 0968
C EXACT CONTAINS THE NG VALUES OF THE SECOND CURVE TO BE PLOTTED THAT 0969
C HAS BEEN COMPUTED IN USERSX (ONLY IF ONLY1=.FALSE.). 0970
C AA IS ( SQUARE ROOT OF THE COVARIANCE MATRIX OF THE SOLUTION)/STDDEV. 0971
C YOU CAN USE AA AND THE LAST 4 ARGUMENTS TO COMPUTE ERROR 0972
C ESTIMATES, AS SHOWN BELOW. 0973
C (CONTIN SETS DOERR=.FALSE. IF IT WAS NOT ABLE TO COMPUTE AA OR 0974
C STDDEV. THE ERROR ESTIMATES THEN CANNOT BE COMPUTED AND ARE SET 0975
C TO ZERO.) 0976
C BELOW IS ILLUSTRATED THE CASE FROM THE ANALYSIS OF CIRCULAR DICHROIC 0977
C SPECTRA, WHERE THE NSPECT-BY-1 VECTOR SOL (THE FRACTIONS OF EACH 0978
C REFERENCE PROTEIN SPECTRUM IN THE SPECTRUM BEING ANALYZED) IS 0979
C RIGHT-MULTIPLIED BY THE NCLASS-BY-NSPECT MATRIX FCAP (THE X-RAY 0980
C VALUES OF THE FRACTIONS OF THE REFERENCE PROTEINS IN EACH CLASS) 0981
C TO YIELD THE NCLASS-BY-1 VECTOR F (THE FRACTION OF THE PROTEIN 0982
C BEING ANALYZED IN EACH CLASS). F AND THE ERROR ESTIMATES ARE 0983
C THEN NORMALIZED BY DIVIDING EACH ELEMENT BY THE SUM OF THE F'S 0984
C (WHICH IS ALSO PRINTED AS THE SCALE FACTOR, SUMF). 0985
C NSPECT = THE NUMBER OF REFERENCE PROTEIN CD SPECTRA. 0986
C NCLASS = THE NUMBER OF CONFORMATIONAL CLASSES. 0987
SUBROUTINE USEROU (CQUAD,G,SOL,EXACT,AA,MG,STDDEV,DOERR,NGLEY) 0988
DOUBLE PRECISION AA 0989
DOUBLE PRECISION PRECIS, RANGE 0990
LOGICAL DOCHOS, DOERR, DOMOM, DOUSIN, DOUSNQ, LAST, 0991
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 0992
DIMENSION CQUAD(MG), G(MG), SOL(MG), EXACT(MG), AA(MG,MG) 0993
C----------------------------------------------------------------------- 0994
C IF YOU CHANGE NSPECT OR NCLASS IN THE DATA STATEMENT BELOW, THEN YOU 0995
C MUST READJUST THE DIMENSIONS IN THE FOLLOWING STATEMENT TO 0996
C DIMENSION FCAP(NCLASS,NSPECT), F(NCLASS), ERROR(NCLASS) 0997
C AND, OF COURSE, MODIFY FCAP (THE X-RAY FRACTIONS) AND NSPECT 0998
C AND NCLASS IN THE DATA STATEMENT BELOW. 0999
C----------------------------------------------------------------------- 1000
DIMENSION FCAP(3,16), F(3), ERROR(3) 1001
COMMON /DBLOCK/ PRECIS, RANGE 1002
COMMON /SBLOCK/ DFMIN, SRMIN, 1003
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 1004
2 SRANGE 1005
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 1006
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 1007
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 1008
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 1009
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 1010
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 1011
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 1012
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 1013
2 LUSER(30) 1014
DATA NSPECT/16/, NCLASS/3/, FCAP/ 1015
1 .79,.00,.21 , .41,.16,.43 , .23,.40,.37 , 1016
2 .28,.14,.58 , .45,.24,.31 , .09,.34,.57 , 1017
3 .02,.51,.47 , .39,.00,.61 , .07,.52,.41 , 1018
4 .24,.15,.61 , .51,.24,.25 , .62,.05,.33 , 1019
5 .37,.15,.48 , .28,.33,.39 , 1020
6 .54,.12,.34 , .26,.44,.30 / 1021
SUMF=0. 1022
DO 110 J=1,NCLASS 1023
F(J)=0. 1024
DO 120 K=1,NSPECT 1025
F(J)=F(J)+SOL(K)*FCAP(J,K) 1026
120 CONTINUE 1027
SUMF=SUMF+F(J) 1028
110 CONTINUE 1029
DO 140 J=1,NCLASS 1030
F(J)=F(J)/SUMF 1031
140 CONTINUE 1032
FACTOR=STDDEV/SUMF 1033
C----------------------------------------------------------------------- 1034
C RUSER(14) IS THE SCALE FACTOR BY WHICH THE DATA SPECTRUM IS 1035
C DIVIDED (SEE USERS MANUAL). NORMALLY, FACTOR=STDDEV ABOVE. 1036
IF (ABS(RUSER(14)) .GT. 0.) SUMF=SUMF/RUSER(14) 1037
C----------------------------------------------------------------------- 1038
C ABOVE, THE NCLASS-BY-1 VECTOR F WAS PRODUCED BY MULTIPLYING THE 1039
C NSPECT-BY-1 SOLUTION SOL BY THE NCLASS-BY-NSPECT MATRIX FCAP. 1040
C THIS IS A VERY COMMON TYPE OF TRANSFORMATION, THAT OCCURS IN 1041
C MANY APPPLICATIONS, E.G., WHEN THE SOLUTION IS REPRESENTED BY 1042
C A SUM OF BASIS (E.G., SPLINE) FUNCTIONS. IN THIS LATTER CASE, 1043
C SOL WOULD BE THE EXPANSION COEFFICIENTS AND F WOULD BE THE 1044
C SOLUTION VALUES AT NCLASS GRID POINTS (SEE USERS MANUAL.) 1045
C THE COMPUTATION OF THE NCLASS-BY-1 VECTOR ERROR, THE 1046
C STANDARD ERROR ESTIMATES FOR F, BELOW IS COMPLETELY GENERAL 1047
C FOR ANY FCAP. SO YOU DO NOT HAVE TO CHANGE ANYTHING BELOW IF 1048
C FCAP, NCLASS, OR NSPECT ARE CHANGED. (USUALLY NSPECT=NGL, BUT 1049
C THIS IS NOT NECESSARY.) 1050
C----------------------------------------------------------------------- 1051
DO 210 J=1,NCLASS 1052
ERROR(J)=0. 1053
IF (.NOT.DOERR) GO TO 210 1054
DO 220 K=1,NGLEY 1055
SUM=0. 1056
DO 230 L=1,NSPECT 1057
SUM=SUM+FCAP(J,L)*AA(L,K) 1058
230 CONTINUE 1059
ERROR(J)=ERROR(J)+SUM**2 1060
220 CONTINUE 1061
ERROR(J)=FACTOR*SQRT(ERROR(J)) 1062
210 CONTINUE 1063
C----------------------------------------------------------------------- 1064
C IF YOU CHANGE THE CLASSES, THEN YOU MUST MODIFY THE FOLLOWING 1065
C FORMAT AND WRITE STATEMENTS. 1066
C----------------------------------------------------------------------- 1067
5140 FORMAT (/23X,5HHELIX,3X,10HBETA-SHEET,4X, 1068
1 9HREMAINDER,16X,12HSCALE FACTOR/ 1069
2 9H FRACTION,F19.2,2F13.2,F28.3/ 1070
3 15H STANDARD ERROR,1P3E13.1) 1071
WRITE (NOUT,5140) F,SUMF,ERROR 1072
RETURN 1073
END 1074
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 1075
C SUBROUTINE USERRG. THIS IS A USER-SUPPLIED ROUTINE (NEEDED 1076
C WHEN NORDER IS INPUT NEGATIVE) TO SET NREG AND TO PUT A 1077
C SPECIAL USER-DEFINED REGULARIZOR IN THE FIRST NREG COLUMNS 1078
C AND ROWS OF REG AND THE RIGHT-HAND-SIDE (R.H.S.) OF THE 1079
C REGULARIZOR IN COLUMN NGLP1 OF REG. 1080
C NOTE - IF IWT = 1 OR 4, THEN USERRG IS CALLED 2 TIMES. 1081
C IF IWT = 2, 3, OR 5, THEN USERRG IS CALLED 4 TIMES. 1082
C THEREFORE, IT IS BEST TO HAVE ANY DATA NEEDED BY USERRG READ 1083
C IN ONLY ONCE AND STORED (E.G., IN RUSER, AS ILLUSTRATED 1084
C BELOW). OTHERWISE, THIS DATA WOULD HAVE TO BE READ 2 OR 4 1085
C TIMES. 1086
C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 1087
C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 1088
C BELOW IS ILLUSTRATED A REGULARIZOR THAT PENALIZES DEVIATIONS OF THE 1089
C SOLUTION FROM AN EXPECTED SOLUTION. THE IDENTITY MATRIX GOES 1090
C INTO THE REGULARIZOR, AND THE EXPECTED SOLUTION IS READ INTO 1091
C RUSER(IUSER(1)),...,RUSER(IUSER(1)+NG-1) AND THEN PUT INTO THE 1092
C R.H.S. OF THE REGULARIZOR (COLUMN NGLP1 OF REG). (SEE 1093
C S. TWOMEY, JACM 10, 97 (1963).) 1094
C LUSER(1) IS INPUT INITIALLY AS .FALSE. AND THEN SET TO .TRUE. 1095
C SO THAT THE DATA IS ONLY READ ONCE. 1096
SUBROUTINE USERRG (REG,MREG,MG,NREG) 1097
DOUBLE PRECISION PRECIS, RANGE 1098
DOUBLE PRECISION REG 1099
LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 1100
1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 1101
DIMENSION REG(MREG,MG) 1102
COMMON /DBLOCK/ PRECIS, RANGE 1103
COMMON /SBLOCK/ DFMIN, SRMIN, 1104
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 1105
2 SRANGE 1106
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 1107
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 1108
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 1109
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 1110
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 1111
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 1112
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 1113
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 1114
2 LUSER(30) 1115
C----------------------------------------------------------------------- 1116
C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 1117
C FOR YOUR APPLICATION. 1118
C----------------------------------------------------------------------- 1119
NREG=NG 1120
DO 110 J=1,NREG 1121
DO 120 K=1,NGL 1122
REG(J,K)=0. 1123
120 CONTINUE 1124
REG(J,J)=1. 1125
110 CONTINUE 1126
J=IUSER(1) 1127
K=J+NG-1 1128
IF (LUSER(1)) GO TO 200 1129
5200 FORMAT (5E15.6) 1130
READ (NIN,5200) (RUSER(L),L=J,K) 1131
WRITE (NOUT,5200) (RUSER(L),L=J,K) 1132
LUSER(1)=.TRUE. 1133
200 IROW=0 1134
DO 210 L=J,K 1135
IROW=IROW+1 1136
REG(IROW,NGLP1)=RUSER(L) 1137
210 CONTINUE 1138
RETURN 1139
END 1140
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 1141
C SUBROUTINE USERSI. THIS IS A USER-SUPPLIED ROUTINE (ONLY 1142
C CALLED WHEN SIMULA IS .TRUE.) FOR CALCULATING EXACT(J) 1143
C (THE SIMULATED DATA BEFORE NOISE IS ADDED) AND 1144
C Y(J) (THE SIMULATED NOISY DATA) FOR J=1,NY. 1145
C EXACT(J) MUST BE COMPUTED IN USEREX. 1146
C IUSER(3) = STARTING INTEGER FOR RANDOM NUMBER GENERATOR RANDOM. 1147
C IUSER(3) AND RUSER(3) MUST BE SUPPLIED BY THE USER. 1148
C IUSER(3) MUST BE BETWEEN 1 AND 2147483646. IF IT IS NOT, THEN 1149
C IT IS SET TO THE DEFAULT VALUE OF 30171. 1150
C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 1151
C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 1152
C BELOW IS ILLUSTRATED THE CASE WHERE ZERO-MEAN 1153
C NORMALLY DISTRIBUTED PSEUDORANDOM NOISE IS ADDED TO EXACT(J). 1154
C SD(J), THE STANDARD DEVIATION OF THE NOISE AT POINT J, IS 1155
C DETERMINED BY IWT AND RUSER(3) AS FOLLOWS - 1156
C IWT = 1 CAUSES SD(J)=RUSER(3) FOR ALL J. 1157
C IWT = 2 CAUSES SD(J)=RUSER(3)*SQRT(EXACT(J)), AS APPROPRIATE 1158
C FOR POISSON STATISTICS. IN THE POISSON CASE, 1159
C RUSER(3) IS JUST A SCALE FACTOR FOR THE CASE THAT 1160
C EXACT(J) IS NOT IN NUMBER OF EVENTS, I.E., 1161
C RUSER(3)=SQRT(EXACT(J)/(NO. OF EVENTS IN CHANNEL J)). 1162
C THUS, EXACT(J)/RUSER(3)**2 IS THE POISSON VARIABLE. 1163
C IF EXACT(J) IS ALREADY IN NO. OF EVENTS, THEN YOU 1164
C SHOULD SET RUSER(3)=1. 1165
C IWT = 3 CAUSES SD(J)=RUSER(3)*EXACT(J). 1166
C IWT = 4 CAUSES SD(J)=RUSER(3)/SQRTW(J). 1167
C IWT = 5 IS FOR THE SPECIAL CASE OF A POISSON SECOND ORDER 1168
C CORRELATION FUNCTION IN PHOTON CORRELATION SPECTROSCOPY. 1169
C YOU SHOULD SET RUSER(3)=1/SQRT(B), AND 1170
C USEREX MUST COMPUTE THE SIMULATED (1ST ORDER CORRELATION 1171
C FUNCTION)*GAMMA, WHERE GAMMA AND B ARE GIVEN IN EQ.(6) 1172
C OF MAKROMOL. CHEM., VOL. 180, P. 201 (1979). 1173
C (SEE ALSO THE USERS MANUAL.) 1174
C----------------------------------------------------------------------- 1175
C CALLS SUBPROGRAMS - USEREX, RGAUSS, ERRMES 1176
C WHICH IN TURN CALL - RANDOM, USERK, USERLF, USERTR 1177
C----------------------------------------------------------------------- 1178
SUBROUTINE USERSI (EXACT,G,MG,MY,SQRTW,T,Y) 1179
DOUBLE PRECISION PRECIS, RANGE 1180
DOUBLE PRECISION DUB 1181
LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 1182
1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 1183
DIMENSION T(MY), EXACT(MY), Y(MY), SQRTW(MY), G(MG) 1184
DIMENSION RN(2), IHOLER(6) 1185
COMMON /DBLOCK/ PRECIS, RANGE 1186
COMMON /SBLOCK/ DFMIN, SRMIN, 1187
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 1188
2 SRANGE 1189
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 1190
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 1191
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 1192
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 1193
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 1194
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 1195
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 1196
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 1197
2 LUSER(30) 1198
DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HS, 1HI/ 1199
TWOPI=6.2831853072D0 1200
DUB=DBLE(FLOAT(IUSER(3))) 1201
IF (DUB.LT.1.D0 .OR. DUB.GT.2147483646.D0) DUB=30171.D0 1202
DO 150 J=1,NY 1203
JJ=J 1204
C----------------------------------------------------------------------- 1205
C RGAUSS DELIVERS TWO NORMAL DEVIATES WITH ZERO MEANS AND STANDARD 1206
C DEVIATIONS 1. THEREFORE IT IS ONLY CALLED FOR ODD J. 1207
C----------------------------------------------------------------------- 1208
K=2-MOD(J,2) 1209
IF (K .EQ. 1) CALL RGAUSS (RN(1),RN(2),TWOPI,DUB) 1210
C----------------------------------------------------------------------- 1211
C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 1212
C FOR YOUR APPLICATION. 1213
C----------------------------------------------------------------------- 1214
EXACT(J)=USEREX(JJ,T,MY,G,MG) 1215
C----------------------------------------------------------------------- 1216
C THE NEXT STATEMENT TEMPORARILY SETS THE ERROR TO A NORMAL DEVIATE 1217
C WITH MEAN = ZERO AND STANDARD DEVIATION = RUSER(3). 1218
C----------------------------------------------------------------------- 1219
ERROR=RUSER(3)*RN(K) 1220
IF (IWT .NE. 2) GO TO 160 1221
IF (EXACT(J) .LT. 0.) CALL ERRMES (1,.TRUE.,IHOLER,NOUT) 1222
ERROR=ERROR*SQRT(EXACT(J)) 1223
GO TO 190 1224
160 IF (IWT .EQ. 3) ERROR=ERROR*EXACT(J) 1225
IF (IWT .NE. 4) GO TO 170 1226
IF (SQRTW(J) .LE. 0.) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) 1227
ERROR=ERROR/SQRTW(J) 1228
170 IF (IWT .NE. 5) GO TO 190 1229
C----------------------------------------------------------------------- 1230
C SPECIAL CASE OF A POISSON SECOND ORDER CORRELATION FUNCTION IN 1231
C PHOTON CORRELATION SPECTROSCOPY. 1232
C G2 = NOISE-FREE NORMALIZED 2ND ORDER CORRELATION FUNCTION. 1233
C G2N1 = NOISY 2ND ORDER CORRELATION FUNCTION LESS 1. 1234
C Y(J) = NOISY FIRST ORDER CORRELATION FUNCTION. IT IS EVALUATED 1235
C BELOW WITH THE FUNCTION SIGN() SO THAT NEGATIVE DATA POINTS 1236
C REMAIN NEGATIVE. THIS PREVENTS A BIAS TOWARD POSITIVE 1237
C NOISE COMPONENTS. 1238
C----------------------------------------------------------------------- 1239
G2=1.+EXACT(J)**2 1240
G2N1=G2-1.+ERROR*SQRT(G2) 1241
Y(J)=SIGN(SQRT(ABS(G2N1)),G2N1) 1242
GO TO 150 1243
190 Y(J)=EXACT(J)+ERROR 1244
150 CONTINUE 1245
RETURN 1246
END 1247
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 1248
C SUBROUTINE USERSX. THIS IS A USER-SUPPLIED ROUTINE (ONLY 1249
C CALLED WHEN YOU HAVE INPUT ONLY1=.FALSE.) 1250
C TO COMPUTE EXACT(J),J=1,NG, WHICH WILL BE PLOTTED WITH 1251
C EACH SOLUTION. USUALLY EXACT IS THE EXACT THEORETICAL 1252
C SOLUTION USED TO SIMULATE YOUR DATA IN USERSI EVALUATED AT 1253
C THE GRID POINTS G(J),J=1,NG. 1254
C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 1255
C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 1256
C BELOW IS ILLUSTRATED THE CASE WHERE EXACT IS 1257
C G(J)**RUSER(8)*EXP(-G(J))/(FACTORIAL OF RUSER(8)), WHERE 1258
C 1. .LE. RUSER(8) .LE. 20. AND THE G(J) ARE NONNEGATIVE. 1259
C----------------------------------------------------------------------- 1260
C CALLS SUBPROGRAMS - GAMLN 1261
C----------------------------------------------------------------------- 1262
SUBROUTINE USERSX (EXACT,G,MG) 1263
DOUBLE PRECISION PRECIS, RANGE 1264
LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 1265
1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 1266
DIMENSION EXACT(MG), G(MG) 1267
COMMON /DBLOCK/ PRECIS, RANGE 1268
COMMON /SBLOCK/ DFMIN, SRMIN, 1269
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 1270
2 SRANGE 1271
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 1272
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 1273
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 1274
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 1275
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 1276
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 1277
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 1278
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 1279
2 LUSER(30) 1280
C----------------------------------------------------------------------- 1281
C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 1282
C FOR YOUR EVALUATION OF EXACT. 1283
C----------------------------------------------------------------------- 1284
IF (RUSER(8).GE.1. .AND. RUSER(8).LE.20.) GO TO 120 1285
5120 FORMAT (/11H RUSER(8) =,E12.4,27H IS OUT OF RANGE IN USERSX.) 1286
WRITE (NOUT,5120) RUSER(8) 1287
STOP 1288
120 EXMIN=-ALOG(SRANGE) 1289
FACTL=GAMLN(RUSER(8)+1.) 1290
DO 150 J=1,NG 1291
EXACT(J)=0. 1292
IF (G(J)) 160,150,180 1293
160 WRITE (NOUT,5160) 1294
5160 FORMAT (/22H NEGATIVE G IN USEREX.) 1295
STOP 1296
180 EX=RUSER(8)*ALOG(G(J))-G(J)-FACTL 1297
IF (EX .GE. EXMIN) EXACT(J)=EXP(EX) 1298
150 CONTINUE 1299
RETURN 1300
END 1301
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 1302
C FUNCTION USERTR. THIS IS A USER-SUPPLIED ROUTINE 1303
C FOR COMPUTING THE QUADRATURE-GRID TRANSFORMATION 1304
C (CALL IT H) H(G) WHEN IFUNCT=1, THE INVERSE TRANSFORMATION 1305
C WHEN IFUNCT=2, AND THE DERIVATIVE OF THE TRANSFORMATION 1306
C WHEN IFUNCT=3. WHEN IGRID=2, G (THE QUADRATURE GRID) WILL 1307
C BE IN EQUAL INTERVALS OF H(G) RATHER THAN IN EQUAL 1308
C INTERVALS OF G (AS IT IS WHEN IGRID=1 AND H IS THE 1309
C IDENTITY TRANSFORMATION). 1310
C BELOW IS ILLUSTRATED THE CASE H(G)=ALOG(G). FOR ANOTHER H, 1311
C YOU CAN REPLACE THE STATEMENTS NUMBERED 210, 220, AND 230. 1312
C THESE ARE THE ONLY STATEMENTS THAT CAN BE REPLACED. ALSO 1313
C NOTE THAT ONLY AN H THAT IS MONOTONIC IN THE RANGE OF 1314
C INTEGRATION MAKES SENSE. 1315
C----------------------------------------------------------------------- 1316
C CALLS SUBPROGRAMS - ERRMES 1317
C----------------------------------------------------------------------- 1318
FUNCTION USERTR (X,IFUNCT) 1319
DOUBLE PRECISION PRECIS, RANGE 1320
LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 1321
1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 1322
COMMON /DBLOCK/ PRECIS, RANGE 1323
COMMON /SBLOCK/ DFMIN, SRMIN, 1324
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 1325
2 SRANGE 1326
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 1327
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 1328
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 1329
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 1330
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 1331
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 1332
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 1333
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 1334
2 LUSER(30) 1335
DIMENSION IHOLER(6) 1336
DATA IHOLER/1HU, 1HS, 1HE, 1HR, 1HT, 1HR/ 1337
IF (IFUNCT.LT.1 .OR. IFUNCT.GT.3) CALL ERRMES (1,.TRUE., 1338
1 IHOLER,NOUT) 1339
IF (IGRID .NE. 1) GO TO 200 1340
C----------------------------------------------------------------------- 1341
C COMPUTE TRANSFORMATION, INVERSE, AND DERIVATIVE FOR IDENTITY 1342
C TRANSFORMATION. 1343
C----------------------------------------------------------------------- 1344
USERTR=1. 1345
IF (IFUNCT .NE. 3) USERTR=X 1346
RETURN 1347
200 IF (IGRID .NE. 2) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) 1348
GO TO (210,220,230),IFUNCT 1349
C----------------------------------------------------------------------- 1350
C COMPUTE TRANSFORMATION. 1351
C----------------------------------------------------------------------- 1352
210 USERTR=ALOG(X) 1353
RETURN 1354
C----------------------------------------------------------------------- 1355
C COMPUTE INVERSE TRANSFORMATION. 1356
C----------------------------------------------------------------------- 1357
220 USERTR=EXP(X) 1358
RETURN 1359
C----------------------------------------------------------------------- 1360
C COMPUTE DERIVATIVE OF TRANSFORMATION. 1361
C----------------------------------------------------------------------- 1362
230 USERTR=1./X 1363
RETURN 1364
END 1365
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 1366
C SUBROUTINE USERWT. THIS IS A USER-SUPPLIED ROUTINE (ONLY 1367
C NEEDED WHEN IWT IS INPUT AS 5) FOR CALCULATING SQRTW (SQUARE 1368
C ROOTS OF THE LEAST SQUARES WEIGHTS) FROM Y, YLYFIT, AND 1369
C ERRFIT, AS EXPLAINED IN DETAIL IN THE USERS MANUAL. 1370
C SEE THE USERS MANUAL FOR THE POSITION IN THE INPUT 1371
C DATA DECK OF ANY INPUT FOR THIS USER SUBPROGRAM. 1372
C BELOW IS ILLUSTRATED THE CASE FROM PHOTON CORRELATION SPECTROSCOPY, 1373
C WHERE THE VARIANCE OF Y IS PROPORTIONAL TO (Y**2+1)/(4*Y**2). 1374
SUBROUTINE USERWT (Y,YLYFIT,MY,ERRFIT,SQRTW) 1375
DOUBLE PRECISION PRECIS, RANGE 1376
LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 1377
1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 1378
DIMENSION Y(MY), YLYFIT(MY), SQRTW(MY) 1379
COMMON /DBLOCK/ PRECIS, RANGE 1380
COMMON /SBLOCK/ DFMIN, SRMIN, 1381
1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 1382
2 SRANGE 1383
COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 1384
1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 1385
2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 1386
3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 1387
4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 1388
5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 1389
COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 1390
1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 1391
2 LUSER(30) 1392
C----------------------------------------------------------------------- 1393
C YOU MUST REPLACE THE FOLLOWING STATEMENTS WITH THOSE APPROPRIATE 1394
C FOR YOUR APPLICATION. 1395
C----------------------------------------------------------------------- 1396
DO 110 J=1,NY 1397
DUM=AMAX1(ABS(Y(J)-YLYFIT(J)),ERRFIT) 1398
SQRTW(J)=2.*DUM/SQRT(DUM*DUM+1.) 1399
110 CONTINUE 1400
RETURN 1401
END 1402
C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 1403
C SUBROUTINE ANALYZ. DOES COMPLETE CONSTRAINED REGULARIZED 1404
C ANALYSIS. 1405
C----------------------------------------------------------------------- 1406
C CALLS SUBPROGRAMS - SETSCA, SEQACC, SETREG, USEREQ, ELIMEQ, SVDRS2, 1407
C ERRMES, DIAREG, DIAGA, SETGA1, SETVAL, LDPETC, RUNRES, 1408
C ANPEAK, SETSGN, SETNNG 1409