SUBROUTINE CGETMOD C C Version 1.0: 29 Dec 1990 J.T. Armstrong C C Adapted for BINFIT, binary data fitting program for Mt. Wilson C interferometer data, from BGETDAT, by Richard S. Simon. C This subroutine interacts with the user to: C 1) Get input file name for existing model, or open a new file C if necessary. C 2) Allows changes to the model to be made. C 3) Allows user to specify which parameters to be varied C during model fit. C C The binary model looks like this: C For each component IC and each filter IG, we have C SEP Separation of component from phase center (mas) C POSANG Position angle East of North (deg) C DELMAG Magnitudes fainter than component 1; for C component 1, zero by default C COLOR Mag(filter IG) minus mag(filter 1) plus color(com- C ponent 1, filter IG) for components IC > 1 C DIAM Diameter (mas) C DRATIO DIAM(IG)/DIAM(1) for modeling limb darkening C FLUX (Relative) flux C FLUXTOT Total flux for filter IG C These are EQUIVALENCEd to the three-dimensional array PARM, for C which the first and second indices are IC and IG, and the third C (usually K in this package) is the parameter number: C K = 1 => Separation C K = 2 => Position angle C K = 3 => Diameter C K = 4 => diameter ratio C K = 5 => Magnitude difference (see below) C K = 6 => Color C K = 7 => Component flux C K = 8 => Total flux (all components) for this filter C C These parameters are further defined or constrained as follows: C Separation: for each component, SEP(IC,IG > 1) = SEP(IC,IG = 1) C i.e., only the separation for IG = 1 is used. C Position angle: as for separation, the value for IG = 1 is C used for all components. C Diameter: For IG = 1, this is entered by the user; for IG > 1, C it is calculated from DIAM(1)*DRATIO(IG). C Diameter ratio: for IG = 1, this is ignored; for IG > 1, it is C DIAM(IC,IG)/DIAM(IC,1) as entered by the user. This allows C the user to fix diameter ratios for different colors according C to some limb-darkening model, and should be set to 1 for no C limb darkening. C Magnitude difference: For IC = 1, this is not used. For C IC > 1, it is -2.5*LOG10( FLUX(IC,IG)/FLUX(1,IG) ), i.e., C the magnitude difference between component IC and component 1 C in this filter. C Color: All colors are relative to filter 1 in the sense C magnitude(filter IG) minus magnitude(filter 1), with colors C for component 1 defined to be zero. Colors for components C IC > 1 are relative to the color of component 1. That is, C if COLOR(IC=1,IG>1) = x and COLOR(IC>1,IG>1) = y, then C DELMAG(IC>1,IG>1) = DELMAG(IC>1,IG=1)+COLOR(IC>1,IG>1)) C so that DELMAG remains the difference in magnitude between C component 1 and component IC in each filter. C Flux: The relative flux of component IC in filter IG. C Total flux: For IC = 1, this is the sum over ALL components IC C of the flux in filter IG. C C The parameters are entered via CREADMDL and are interpreted as follows: C For any line, first number is component, second is wavelength (nm) C IC = 1 and IG = 1: C Third number is separation; C Fourth is position angle; C Fifth is diameter; C Sixth is total flux in filter 1. C NOTE: magnitude difference and flux of component 1 cannot C be entered. C IC = 1 and IG > 1: C Third number is diameter ratio. C Fourth number is total flux in filter IG. C NOTE: no other parameters may be entered. C IC > 1 and IG = 1: C Third number is separation; C Fourth is position angle; C Fifth is diameter; C Sixth is magnitude difference from component 1, filter 1. C NOTE: magnitude difference and total flux in filter 1 cannot C be entered. C IC > 1 and IG > 1: C Third number is the diameter ratio; C Fourth number is color in filter IG relative to filter 1. C NOTE: no other parameters may be entered. C C C The variation of various parameters is determined by their inter- C relationships. The user will expect that telling the program to vary C d(mag) for one parameter will reduce the number of degrees of freedom by C one. This is implemented by fixing the total flux; thus varying d(mag) C for filter 1 for one component will cause all the fluxes of that C component and of component 1 to change. However, if the user specifies C that the total flux and d(mag) for one component should change s/he will C expect that the degrees of freedom will be reduced by two. This is C implemented by allowing the fluxes of one component to vary C independently from those of the other, while fluxes within one component C stay at fixed ratios. If a color is allowed to vary as well, another C constraint is lifted by allowing flux ratios within a component to vary. C A similar situation will apply to the diameters; however, given the Mark C III's limitations, it is probably unsafe to try to determine the C (relative) amount of limb darkening in the different channels; it is C best just to put in the appropriate values from your favorite model. C C C Thus, when the user is asked to indicate which parameters are to C vary, some care is needed in determining which ones actually will C change, and how. Here, "vary" means that the parameter is one whose C best value is sought; "change" means that the value of the parameter C may change because some other parameter has changed. C C IC = 1 and IG = 1: C Separation and position angle may each vary, but will C normally both be zero; C Diameter may vary; C Diameter ratio may NOT vary (it isn't used); C Magnitude difference may NOT vary (it isn't used); C Color may NOT vary (and the color is always zero); C Flux may NOT vary; it may change if the total flux and C the magnitude difference for another component vary; C Total flux may vary; if it does: C If no magnitude differences for other components vary, C then fluxes for all components change together; C If some magnitude difference or color for another component C varies, then the flux for that component and the C total flux will vary independently, which means C that the flux for component 1 will change. C C IC = 1 and IG > 1: C Separation and position angle may NOT vary; C Diameter may NOT vary, but it may change if the diameter C ratio varies; C Diameter ratio DIAM(1,IG)/DIAM(1,1) may vary, but it's C not a good idea (see above); C Magnitude difference may NOT vary; C Color may NOT vary (and the color is always zero); C Flux and total flux may change or vary under the conditions C given for IG = 1. C C IC > 1 and IG = 1: C Separation and position angle may each vary; C Diameter may vary; C Diameter ratio may NOT vary (it isn't used); C Magnitude difference may vary: C If the total flux for this filter does not vary, the C fluxes of this component and of component 1 will C change together; C If the total flux does vary, the flux of this compon- C ent and the total flux will vary independently, C which means that the flux for component 1 will change; C Color may NOT vary (it isn't used); C Flux may NOT vary; it will change if the magnitude difference C varies; C Total flux may NOT vary (it isn't used). C C IC > 1 and IG > 1: C Separation and position angle may NOT vary; C Diameter may NOT vary, but it may change if the diameter C ratio varies; C Diameter ratio DIAM(IC,IG)/DIAM(IC,1) may vary, but it's C not a good idea; C Magnitude difference may NOT vary, but it may change if the C color varies; C Color may vary: C If the total flux for this filter does not vary, the C flux of this component and of component 1 will C change together; C If the total flux does vary, the flux of this compon- C ent and the total flux will vary independently, C which means that the flux for component 1 will change; C Flux may NOT vary; it will change if the magnitude difference C varies; C Total flux may NOT vary (it isn't used). C----------------------------------------------------------------------- INCLUDE 'BINFIT.INC' INCLUDE 'MODEL.INC' CALL PUTOUT(' ') CALL GETMNAME ! Get model file name D IF (NEWMOD) THEN D WRITE(OUTC,*) ' NEWMOD is true.' D ELSE IF(.NOT.NEWMOD) THEN D WRITE(OUTC,*) ' NEWMOD is false.' D ELSE D WRITE(OUTC,*) ' NEWMOD not set true or false.' D END IF IF (NEWMOD) THEN CALL MNEWMOD ! Make new model file ELSE CALL ROLDMOD ! Read pre-existing model file END IF CALL CADJMOD ! Adust model fluxes C 65 WRITE (OUTC,1060) ! Return here after changes CALL CWRITEMD(OUTC,' ',PARM,0,0) CALL PUTOUT(' ') IF (.NOT.QUERY(' Is this acceptable? (Y or N): ')) THEN WRITE (OUTC,1080) READ (INC, 1000) ANSWER CALL UPCASE (ANSWER) IF (ANSWER .EQ. 'C' .OR. ANSWER .EQ. 'B') THEN CALL CHNGMOD ! Changes to model CALL CADJMOD ! Recalculate fluxes END IF IF (ANSWER .EQ. 'A' .OR. ANSWER .EQ. 'B') THEN CALL ADDMOD ! Additions to model CALL CADJMOD ! Recalculate fluxes END IF ! End Additions/Changes CLOSE (UNIT=INMOD, STATUS='delete') OPEN (UNIT=INMOD, FILE=MODDSN, STATUS='NEW',FORM='FORMATTED') CALL CWRITEMD(OUTC,' ',PARM,0,0) C Zero all other model components ISTART = NCOMP + 1 DO IC = ISTART, MXMOD DO IG = 1, NFILT DELMAG(IC,IG) = DMAGMAX FLUX(IC,IG) = 0. END DO END DO GO TO 65 ! Return to ask if there are more changes END IF NPARM = NCOMP*NFILT*8 D WRITE(OUTC,'(A,I3)') ' In CGETMOD: NPARM = ',NPARM C Format statements 1000 FORMAT (A,A,A,A) 1060 FORMAT (//' Starting Model:') 1080 FORMAT ('$Do you wish to make Changes/Deletions, Additions,', + ' or Both (C, A, or B) ? ') 1090 FORMAT (' Changes to ',I2,')',4X,7F10.4) 1160 FORMAT (' Enter bandwidth for ',F3.0,' nm filter.') C Ruler line for output: C!II.LLL..+SS.SSSf.+PAP.PAPf.+D.DIf.+.DRDf.(+D.DMAf.+C.COLf.(+F.FLU).+F.TOTf write(outc,*)' You can now enter radius and position angle ', 1'gradients for comp. 2' write(outc,*)' in mas/hour and degree/hour!' write(outc,*)' Enter 0 0, if you do not want to use this option.' read(inc,*)dsep,dposang c Convert dsep, dposang to rad/day dsep=dsep*24.*mas2rad dposang=dposang*24.*degrad RETURN END