C++++++++++++++++ DOUBLE PRECISION VERSION 2DP (MAR 1984) ++++++++++++++ 0003 C CONTIN. MAIN SUBPROGRAM. 0004 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 AS A NEW USER YOU SHOULD NOTIFY S.W. Provencher at: C sp@S-provencher.COM C C SO THAT YOU CAN OBTAIN THE USERS MANUAL (WHICH IS ESSENTIAL) 0019 C AND HAVE YOUR NAME PUT ON A MAILING LIST FOR UPDATES, ETC. 0020 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. 50) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) 0708 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 C WHICH IN TURN CALL - H12, GETROW, USERK, USERLF, USERRG, LH1405, 1410 C QRBD, G1, G2, DIFF, LDP, CVNEQ, FISHNI, BETAIN, GAMLN, 1411 C NNLS, PLRES, UPDSGN, UPDDON, FFLAT, UPDLLS, 1412 C GETPRU, GETYLY, PGAUSS, MOMENT, MOMOUT, PLPRIN, USEROU, USERTR 1413 C----------------------------------------------------------------------- 1414 SUBROUTINE ANALYZ (ISTAGE, 1415 1 A,AA,AEQ,AINEQ,CQUAD,EXACT,G,IISIGN,IWORK,LBIND,LSDONE,MA, 1416 2 MDONE,MEQ,MG,MINEQ,MREG,MWORK,MY,PIVOT,REG,RHSNEQ,S,SOLBES, 1417 3 SOLUTN,SQRTW,SSCALE,T,VALPCV,VALPHA,VDONE,VK1Y1,WORK, 1418 4 Y,YLYFIT) 1419 DOUBLE PRECISION PRECIS, RANGE 1420 DOUBLE PRECISION A, AA, ABS, AEQ, AINEQ, ALPBES, ALPHA, 1421 1 ALPOLD, DUB, ONE, PIVOT, RALPFL, REG, RHSNEQ, 1422 2 S, SOLBES, SOLUTN, SSCALE, VALPCV, VALPHA, VK1Y1, WORK, ZERO 1423 LOGICAL DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, NEWPG1, 1424 1 NONNEG, ONLY1, PRWT, PRY, SIMULA, LUSER 1425 LOGICAL LDUM, LBIND, FLAT, PPLTPR, HEADNG, NEWPAG, NOPAGE 1426 DIMENSION A(MA,MG), T(MY), Y(MY), SQRTW(MY), G(MG), CQUAD(MG), 1427 1 REG(MREG,MG), AEQ(MEQ,MG), PIVOT(MEQ), VK1Y1(MG), S(MG,3), 1428 2 AINEQ(MINEQ,MG), VALPHA(MG), VALPCV(MG), RHSNEQ(MINEQ), 1429 3 WORK(MWORK), IWORK(MA), EXACT(MG), SOLUTN(MG), LBIND(MINEQ), 1430 4 IISIGN(MG), SOLBES(MG), LSDONE(MDONE,3,2), VDONE(MDONE), 1431 5 YLYFIT(MY), AA(MG,MG), SSCALE(MG) 1432 DIMENSION IHOLER(6), PREJ(2), LLSIGN(5), RS2MNX(2) 1433 COMMON /DBLOCK/ PRECIS, RANGE 1434 COMMON /SBLOCK/ DFMIN, SRMIN, 1435 1 ALPST(2), EXMAX, GMNMX(2), PLEVEL(2,2), RSVMNX(2,2), RUSER(551), 1436 2 SRANGE 1437 COMMON /IBLOCK/ IGRID, IQUAD, IUNIT, IWT, LINEPG, 1438 1 MIOERR, MPKMOM, MQPITR, NEQ, NERFIT, NG, NINTT, NLINF, NORDER, 1439 2 IAPACK(6), ICRIT(2), IFORMT(70), IFORMW(70), IFORMY(70), 1440 3 IPLFIT(2), IPLRES(2), IPRINT(2), ITITLE(80), IUSER(50), 1441 4 IUSROU(2), LSIGN(4,4), MOMNMX(2), NENDZ(2), NFLAT(4,2), NGL, 1442 5 NGLP1, NIN, NINEQ, NNSGN(2), NOUT, NQPROG(2), NSGN(4), NY 1443 COMMON /LBLOCK/ DOCHOS, DOMOM, DOUSIN, DOUSNQ, LAST, 1444 1 NEWPG1, NONNEG, ONLY1, PRWT, PRY, SIMULA, 1445 2 LUSER(30) 1446 DATA IHOLER/1HA, 1HN, 1HA, 1HL, 1HY, 1HZ/, RALPHA/1./ 1447 ABS(ALPHA)=DABS(ALPHA) 1448 SINGLE(ONE)=SNGL(ONE) 1449 C SINGLE(VAR)=VAR!SP 1450 C ZERO=0.E0!SP 1451 ZERO=0.D0 1452 C ONE=1.E0!SP 1453 ONE=1.D0 1454 NGLE=NGL-NEQ 1455 C----------------------------------------------------------------------- 1456 C PUT SCALE FACTORS FOR SOLUTION IN SSCALE AND SCALE INEQUALITY 1457 C CONSTRAINTS. 1458 C----------------------------------------------------------------------- 1459 MG1=MG+1 1460 MG2=MG1+MG 1461 MG3=MG2+MG 1462 CALL SETSCA (WORK,WORK(MG1),WORK(MG2),WORK(MG3), 1463 1 AINEQ,CQUAD,G,ISTAGE,MG,MINEQ,MREG,MY,NGLE,REG,S,SQRTW,SSCALE,T, 1464 2 Y) 1465 C----------------------------------------------------------------------- 1466 C DO SEQUENTIAL HOUSEHOLDER TRANSFORMATIONS TO COMPRESS 1467 C COEFFICIENT MATRIX INTO AN UPPER TRIANGLE. 1468 C----------------------------------------------------------------------- 1469 CALL SEQACC ( 1470 1 A,CQUAD,G,IUNIT,IWT,MA,MG,NG,NGL,NGLP1,NLINF,NY, 1471 2 RANGE,SQRTW,SSCALE,T,Y) 1472 NGLY=MIN0(NGL,NY) 1473 C----------------------------------------------------------------------- 1474 C SET UP REGULARIZOR. 1475 C----------------------------------------------------------------------- 1476 CALL SETREG (MG,MREG,NENDZ,NG,NGL,NGLE,NGLP1,NORDER, 1477 1 NOUT,NREG,REG,SSCALE) 1478 IF (NEQ .LE. 0) GO TO 200 1479 C----------------------------------------------------------------------- 1480 C SET EQUALITY CONSTRAINTS IN USER-SUPPLIED PROGRAM USEREQ AND THEN 1481 C SCALE THEM AND NORMALIZE THEIR ROWS (IN THE L1 METRIC). 1482 C----------------------------------------------------------------------- 1483 CALL USEREQ (AEQ,CQUAD,MEQ,MG) 1484 DO 150 ICOL=1,NGL 1485 DO 160 IROW=1,NEQ 1486 AEQ(IROW,ICOL)=AEQ(IROW,ICOL)*SSCALE(ICOL) 1487 160 CONTINUE 1488 150 CONTINUE 1489 DO 170 IROW=1,NEQ 1490 DUB=ZERO 1491 DO 180 ICOL=1,NGL 1492 DUB=DUB+ABS(AEQ(IROW,ICOL)) 1493 180 CONTINUE 1494 IF (DUB .LE. ZERO) CALL ERRMES (0,.TRUE.,IHOLER,NOUT) 1495 DUB=ONE/DUB 1496 DO 190 ICOL=1,NGLP1 1497 AEQ(IROW,ICOL)=AEQ(IROW,ICOL)*DUB 1498 190 CONTINUE 1499 170 CONTINUE 1500 C----------------------------------------------------------------------- 1501 C ELIMINATE EQUALITY CONSTRAINTS USING A BASIS OF NULL SPACE OF AEQ. 1502 C----------------------------------------------------------------------- 1503 200 CALL ELIMEQ (AEQ,MEQ,MG,PIVOT,NEQ,NGL,A,MA,REG,MREG,NREG,NGLP1, 1504 1 NGLY,VK1Y1,RANGE) 1505 C----------------------------------------------------------------------- 1506 C DO SINGULAR VALUE DECOMPOSITION AND DIAGONALIZATION OF REGULARIZOR. 1507 C----------------------------------------------------------------------- 1508 CALL SVDRS2 (REG(1,NEQ+1),MREG,NREG,NGLE,REG(1,NGLP1),MREG, 1509 1 1,S,IERROR,RANGE) 1510 IF (IERROR .NE. 1) CALL ERRMES (1,.TRUE.,IHOLER,NOUT) 1511 CALL DIAREG (A,AEQ,MA,MEQ,MG,MREG,NEQ,NGL,NGLE,NGLY,NOUT, 1512 1 NUNREG,PIVOT,PRECIS,RANGE,REG,S) 1513 C----------------------------------------------------------------------- 1514 C DO SINGULAR VALUE DECOMPOSITION AND DIAGONALIZATION OF MODIFIED 1515 C COMPRESSED COEFFICIENT MATRIX A. 1516 C----------------------------------------------------------------------- 1517 CALL SVDRS2 (A,MA,NGLY,NGLE,A(1,NGLP1),MA,1,S,IERROR,RANGE) 1518 IF (IERROR .NE. 1) CALL ERRMES (2,.TRUE.,IHOLER,NOUT) 1519 CALL DIAGA (A,MA,MG,MREG,NEQ,NGL,NGLE,NGLP1,NGLY,REG,S) 1520 DO 320 J=1,NGLE 1521 YLYFIT(J)=S(J,1) 1522 320 CONTINUE 1523 5320 FORMAT (/16H SINGULAR VALUES/(1X,1P10E13.3)) 1524 WRITE (NOUT,5320) (YLYFIT(J),J=1,NGLE) 1525 C----------------------------------------------------------------------- 1526 C LEFT MULTIPLY BY INEQUALITY MATRIX. 1527 C----------------------------------------------------------------------- 1528 IF (NINEQ .GT. 0) CALL SETGA1 (NINEQ, 1529 1 A,AINEQ,MA,MG,MINEQ,MREG,NGL,NGLE,REG) 1530 IF (MAX0(NQPROG(1),NQPROG(2)).LE.0 .AND. ALPST(ISTAGE).LE.0.) CALL 1531 1 ERRMES (3,.TRUE.,IHOLER,NOUT) 1532 RSVM2J=SRANGE 1533 IF (DFMIN.LT.0. .OR. DFMIN+FLOAT(NUNREG).GE.FLOAT(NGLE)) GO TO 390 1534 C----------------------------------------------------------------------- 1535 C GUARANTEE AT LEAST NUNREG+DFMIN DEGREES OF FREEDOM (IGNORING 1536 C INEQUALITY CONSTRAINTS) BY USING RSVM2J INSTEAD OF RSVMNX(2,J), 1537 C IF NECESSARY. 1538 C----------------------------------------------------------------------- 1539 ALP2=S(1,1)**2 1540 DO 350 J=1,500 1541 IF (ALP2 .LE. 0.) GO TO 390 1542 DF=0. 1543 DO 360 K=1,NGLE 1544 DUM=S(K,1)**2 1545 DF=DF+DUM/(DUM+ALP2) 1546 360 CONTINUE 1547 IF (DF .GT. DFMIN+FLOAT(NUNREG)) GO TO 370 1548 ALP2=ALP2*.1 1549 350 CONTINUE 1550 370 RSVM2J=SQRT(10.*ALP2)/S(1,1) 1551 390 BTEST=SRANGE 1552 ALPBES=ZERO 1553 VARZ=SRANGE 1554 RS2MNX(1)=-1. 1555 RS2MNX(2)=-1. 1556 C----------------------------------------------------------------------- 1557 C START OF NQPROG(1) REGULARIZED SOLUTIONS 1558 C----------------------------------------------------------------------- 1559 LDUM=.TRUE. 1560 PPLTPR=MAX0(IPRINT(ISTAGE),IUSROU(ISTAGE)) .GE. 3 1561 NOPAGE=.NOT.PPLTPR 1562 LHEDNG=0 1563 IF (NOPAGE) LHEDNG=4 1564 K=NQPROG(1) 1565 IF (K .LE. 1) GO TO 410 1566 RTOT=AMIN1(RSVM2J,RSVMNX(2,1))/(RSVMNX(1,1)*PRECIS) 1567 IF (RTOT .LE. 1.) RTOT=RSVMNX(2,1)/(RSVMNX(1,1)*PRECIS) 1568 IF (RTOT .LE. 1.) CALL ERRMES (4,.TRUE.,IHOLER,NOUT) 1569 RALPHA=RTOT**(1./FLOAT(K-1)) 1570 410 ALPHA=RSVMNX(1,1)*PRECIS*S(1,1) 1571 DO 420 J=1,K 1572 CALL SETVAL (ALPHA,LDUM,NINEQ, 1573 1 A,AINEQ,MA,MG,MINEQ,MREG,NGL,NGLE,REG,RHSNEQ,S,VALPCV,VALPHA, 1574 2 VK1Y1) 1575 NEWPAG=MAX0(IPRINT(ISTAGE),IUSROU(ISTAGE), 1576 1 IPLRES(ISTAGE)+1).GE.4 .OR. LDUM 1577 HEADNG=PPLTPR .OR. LDUM .OR. 1578 1 MAX0(IPLRES(ISTAGE),IPLFIT(ISTAGE)).GE.3 1579 CALL LDPETC (3,.TRUE.,NINEQ,.TRUE.,ICRIT(ISTAGE),DOMOM,PPLTPR, 1580 1 .TRUE.,ALPHA,HEADNG,NEWPAG,ALPBES,VAR, 1581 2 A,AA,AINEQ,BTEST,CQUAD,DEGFRE,DEGFRZ,EXACT,G,IERROR, 1582 3 ISTAGE,IWORK,LBIND,MA,MG,MINEQ,MREG,MWORK,MY, 1583 4 NGLE,NGLY,PREJ,REG,RHSNEQ,RS2MNX,S,SOLBES, 1584 5 SOLUTN,SQRTW