C A CODE FOR Teff-LOG G DETERMINATION C FROM H5 AND H6 BALMER LINE PROFILES C VERSION OCT 15., 1996 C c authors: J. Budaj c contact: budaj@ta3.sk c published: http://www.ta3.sk/~budaj/software c special requirements: DBOS environement c non original routins: intep, inte2 c (Press.et.al. 1992, Numerical recipes) c input: keyboard + file c output: balmfit.out c C THE CODE IS WRITTEN FOR THE UNIV. OF SALFORD COMPILER OPERATING C UNDER THE DOS JUST TO ALLOWE YOU TO USE THE GRAPHICS C AND SEE THE FIT. HOWEVER, THE GRAFICS ROUTINES CAN BE SIMPLY C REMOVED AND THEY ARE MARKED IN THE SOURCE FILE THROUGH "!!" C OPTIONS(SAVE,DREAL) C OPTIONS(SAVE,DREAL,CHECK) C !!REMOVE THIS LINE IF YOU USE OTHER COMPILER OPTIONS(SAVE,DREAL,OPTIMISE) C C NMOM1 - NUMBER OF MODELS(+20 INTERP.MODELS), C NMOM - REAL NUMBER OF MODELS(+20 INTERP.MODELS), C NMOD - NUMBER OF GRID MODELS C NSYN - NUMBER OF TABULATED POINTS IN BALM. PROFIL C NSYN2 = 2*NSYN-1 !!! IT MUST BE !!! C NOBS - NUMBER OF OBSERVED POINTS IN BALM. PROFIL C NTEFF - NUMBER OF POINTS IN TEMPERATURE GRID C NLOGG - NUMBER OF POINTS IN LOG G GRID C NEZ2 - NUMBER OF INTERPOL. MODELS PARAMETER (NMOM1=103,NSYN=26,NSYN2=51,NOBS=5000, P NTEFF=18,NLOGG=4,NEZ2=20) IMPLICIT DOUBLE PRECISION(A-H,O-Z) CHARACTER*12 INN1,INN2 DIMENSION CM(NOBS),B(NOBS),TEP(NMOM1),GRAV(NMOM1) P ,H(NMOM1,NSYN),CC(NSYN),TINTRP(NEZ2),GINTRP(NEZ2) P ,OMM(NMOM1),LL(NMOM1),D(NMOM1) P ,GRIDT(NTEFF),GRIDG(NLOGG),PROF(NTEFF,NLOGG) P ,HMIRR(NSYN2),CMIRR(NSYN2) C INPUT DIALOG WRITE(*,*)' BEFORE YOU RUN THIS CODE:' WRITE(*,*)' DOS-COMMAND ' WRITE(*,*)' ALLOWS YOU TO PRINT FIGURES' WRITE(*,*)' BY MEANS OF [PRTSC]' WRITE(*,*)' EXAMPLE OF INPUT FILE: BALMFIT.IN' OPEN(9,FILE='BALMFIT.OUT',STATUS='UNKNOWN') 10 WRITE(*,*) WRITE(*,'(A,$)')' INPUT DATA FILE => ' READ(5,'(A12)')INN1 OPEN(7,FILE=INN1,STATUS='OLD') C WRITE(*,*) C WRITE(*,*)'---------------------' C WRITE(*,*)' START -1- ' C WRITE(*,*)' END -2- ' C WRITE(*,'(A,$)')' --------------------- => ' C READ(5,*)IAD C IF (IAD.EQ.2)GOTO 400 C------------------------------SPECTROSCOPY----------------------------- C READS FROM THE INPUT FILE READ(7,*) READ(7,*) I=1 READ(7,*)IEPER,CENTER,UNIT READ(7,*) 20 READ(7,*,END=30)CM(I),B(I) CM(I)=(CM(I)-CENTER)*UNIT B(I)=B(I)*1.0D3 I=I+1 GOTO 20 30 NBOD=I-1 WRITE(*,*) WRITE(*,'(30X,A15,I5)')' No. OF POINTS ',NBOD CLOSE(7) WRITE(*,*) WRITE(*,*)'SCATTERED LIGHT IN %,' WRITE(*,'(A,$)')' IF NO PUT -0.- => ' READ(5,*)AM C OBSERVED PROFIL WILL BE CORRECTED TO THE SCATTERED LIGHT IF (AM.GT.0.)THEN DO 60 I=1,NBOD B(I)=(B(I)-AM/100.*1000.)/(1.-AM/100.) 60 CONTINUE ENDIF 15 WRITE(*,*) WRITE(*,*)'TO NORMALIZE SYNTHETIC KURUCZ PROFILES GIVE AT WHICH' WRITE(*,*)'DELTA LAMBDA [NM] THEY SHOULD MERGE INTO COMTINUUM' WRITE(*,'(A,$)')' IF NO PUT -0.- => ' READ(5,*)AL WRITE(*,*) WRITE(*,*)'ROTATION V*SIN I [KM/S]' WRITE(*,'(A,$)')' IF NO PUT -0.- => ' READ(5,*)BV IF (BV.GT.0.) THEN WRITE(*,*) WRITE(*,*)'DOUBLE LETTER CODE OF EXPECTED SP. TYPE "XY"' WRITE(*,*)'(X: O=0,B=1,A=2,F=3,G=4,K=5, Y=No. OF SP.SUBCLASS' WRITE(*,'(A,$)')' E.G. A4=24) => ' READ(5,*)IBA ENDIF WRITE(*,*) WRITE(*,*)'SOLAR ABUNDANCES -0- ' WRITE(*,*)'METALS 10-TIMES LOWERED -1- ' WRITE(*,'(A,$)')' METALS 100-TIMES LOWERED -2- => ' READ(5,*)IEPP WRITE(*,*) WRITE(*,*)'TEMPERATURE RANGE OF INTEREST T1 => ' READ(5,*)CA,CB,CD,CE WRITE(*,*) WRITE(*,*)'HOW MANY ADDITIONAL NON-GRID MODELS SHALL I CONSIDER' WRITE(*,'(A,$)')' AND INTERPOLATE [0 - NONE, 20 -MAXIM.] => ' READ(5,*)INTRP IF (INTRP.GE.1)THEN DO 17 I=1,INTRP WRITE(*,*)' RANGE ALLOWED: 7000 OF INTERP. MODELS => ' READ(5,*)TINTRP(I),GINTRP(I) 17 CONTINUE ENDIF WRITE(*,*) WRITE(*,'(15X,A,A)')' NOTE: OBSERVED PROFIL IS CORRECTED FOR', P ' SCATTERED LIGHT' WRITE(*,'(15X,A)')' ALL SYNTHETIC PROFILES ARE NORMALIZED' WRITE(*,'(15X,A,A)')' ONLY SYNTH. PROF. WITHIN GIVEN', P ' INTERVAL ARE CORRECTED' WRITE(*,'(15X,A)')' FOR ROTATION' C WRITE(*,*)' TABULATED DELTA LAMBDA' C WRITE(*,*)' 0,0.02,0.04,0.06,0.08,0.1,0.15,0.2,' C WRITE(*,*)' 0.3,0.4,0.5,0.6,0.8,1,1.2,1.4,1.6,' C WRITE(*,*)' 1.8,2,2.2,2.4,2.8,3.2,3.6,4,5 NM' C WRITE(*,*) WRITE(*,*) WRITE(*,*)'---------------------' WRITE(*,*)' REPEATE DIALOGUE FROM BEGINNING -1- ' WRITE(*,*)' REPEATE IT WITHOUT LOADING DATA -2- ' WRITE(*,*)' END -3- ' WRITE(*,*)' GO ON -4- ' WRITE(*,'(A,$)')' --------------------- => ' READ(5,*)IAD IF (IAD.EQ.1)GOTO 10 IF (IAD.EQ.2)GOTO 15 IF (IAD.EQ.3)GOTO 400 C---------------THIS READS SYNTHETIC PROFILES--------------------------- DATA CC/0.,0.02,0.04,0.06,0.08,0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.8 P ,1.,1.2,1.4,1.6,1.8,2.,2.2,2.4,2.8,3.2,3.6,4.,5./ CALL SYNPRO(TEP,GRAV,H,NMOM1,NMOM,NSYN,INTRP,NMOD,IEPER,IEPP) C------------INTERPOLATION OF ADDITIONAL SYNTHETIC PROFILES------------- DO 100 I=1,INTRP TEP(NMOD+I)=TINTRP(I) GRAV(NMOD+I)=GINTRP(I) 100 CONTINUE DATA GRIDT/7000.,7500.,8000.,8500.,9000.,9500.,10000.,11000., P 12000.,13000.,14000.,15000.,16000.,17000.,18000.,20000., P 22500.,25000./ DATA GRIDG/3.,3.5,4.,4.5/ IF (INTRP.GE.1)THEN DO 170 I=1,INTRP DO 160 J=1,NSYN DO 150 K=1,NTEFF DO 140 L=1,NLOGG PROF(K,L)=0.0 DO 130 M=1,NMOD IF(DABS(TEP(M)-GRIDT(K)).LT.1.D-30.AND. P DABS(GRAV(M)-GRIDG(L)).LT.1.D-30)THEN PROF(K,L)=H(M,J) GOTO 130 ENDIF 130 CONTINUE IF(DABS(PROF(K,L)-0.0).LT.1.D-30)THEN WRITE(*,*)' T-LOGG GRID IS NOT COMPLETE TO INTERP.' PAUSE ENDIF 140 CONTINUE 150 CONTINUE CALL INTEP2(GRIDT,GRIDG,PROF,NTEFF,NLOGG,TINTRP(I), P GINTRP(I),H(NMOD+I,J)) 160 CONTINUE 170 CONTINUE ENDIF C---------------ROTATIONAL BROADENING OF MODELS WITHIN GIVEN INTERVAL--- IF (BV.GT.0.)THEN C CALCULATES DOPPLER WIDTH -DLL, C LIMB. DARK. COEFF. -UKOEF CALL STEM(IBA,IEPER,BV,UKOEF,DLL) C CREATES EVEN NUMBER OF STEPS (NSPACE) IN FILTER WIDTH SPACE0=DLL/20. IF (SPACE0.GT.CC(2)-CC(1))SPACE0=CC(2)-CC(1) NSPACE=DINT(2.D0*DLL/SPACE0)+1 IF (DBLE(NSPACE)/2.D0-DBLE(NSPACE/2).GT.0.1) NSPACE=NSPACE+1 SPACE=2.D0*DLL/DBLE(NSPACE) C WRITE(5,'(/15X,A,A,A,E11.3,I5,E11.3)')'FILTER WIDTH, ', C P ' No. OF STEPS IN IT, ',' MODIF.STEP', C P 2.D0*DLL,NSPACE,SPACE DO 80 I=1,NMOM IF (TEP(I).GE.CA.AND.TEP(I).LE.CB.AND.GRAV(I).GE.CD.AND. P GRAV(I).LE.CE)THEN HMIRR(NSYN)=H(I,1) CMIRR(NSYN)=CC(1) DO 78 J=2,NSYN HMIRR(J-1)=H(I,NSYN-J+2) CMIRR(J-1)=-CC(NSYN-J+2) HMIRR(J-1+NSYN)=H(I,J) CMIRR(J-1+NSYN)=CC(J) 78 CONTINUE WRITE(*,'(F6.0,F6.2,A)')TEP(I),GRAV(I), P ' IS BEING CORRECTED FOR ROTATION' DO 79 J=1,NSYN H(I,J)=AKONV(CC(J),CMIRR,HMIRR,SPACE,NSPACE, P DLL,NSYN2) 79 CONTINUE ENDIF 80 CONTINUE ENDIF C---------------NORMALIZATION OF ALL SYNTH. PROFILES-------------------- IF (AL.GT.0.)THEN DO 90 I=1,NMOM CALL NORM(I,NMOM1,NSYN,AL,CC,H) 90 CONTINUE ENDIF C---------------COMPARING OF OBSERV. AND SYNTH. PROFILES---------------- IKK=0 DO 210 I=1,NMOM D(I)=0. 210 CONTINUE DO 250 I=1,NMOM IF (TEP(I).GE.CA.AND.TEP(I).LE.CB.AND.GRAV(I).GE.CD.AND. P GRAV(I).LE.CE) THEN IKK=IKK+1 IKK1=0 DO 230 IL1=1,NSYN KRA=NSYN-IL1+1 IF(-CC(KRA).GT.CM(1).AND.-CC(KRA).LT.CM(NBOD))THEN CALL INTEP(-CC(KRA),YDEL,CM,B,NBOD,IER) D(I)=(H(I,KRA)-YDEL)*(H(I,KRA)-YDEL)+D(I) IKK1=IKK1+1 ENDIF IF(CC(IL1).GT.CM(1).AND.CC(IL1).LT.CM(NBOD))THEN CALL INTEP(CC(IL1),YDEL,CM,B,NBOD,IER) D(I)=(H(I,IL1)-YDEL)*(H(I,IL1)-YDEL)+D(I) IKK1=IKK1+1 ENDIF 230 CONTINUE D(I)=SQRT(D(I))/DBLE(IKK1) ENDIF 250 CONTINUE C----------------------------------------------------------------------- WRITE(*,*) WRITE(*,'(A,$)')' HOW MANY BEST MODELS AM I CHOOSE => ' READ(5,*)IVOL WRITE(*,*) IF (IKK.LE.IVOL) IVOL=IKK CALL ZORAD(D,OMM,LL,NMOM,IVOL,NMOM1) WRITE(*,*) WRITE(9,*) WRITE(9,'(15X,A21,I1)')' FROM BALMER LINE H',IEPER WRITE(9,*) DO 300 J=1,IVOL WRITE(*,310)LL(J),TEP(LL(J)),GRAV(LL(J)),OMM(J) WRITE(9,310)LL(J),TEP(LL(J)),GRAV(LL(J)),OMM(J) 300 CONTINUE 310 FORMAT(' No:',I2,' T=',F6.0,' LOG G=',F4.2, P ' DEV.:',E11.4) 320 WRITE(*,*) WRITE(*,*)'------------------' WRITE(*,*)'VIEW -1- ' WRITE(*,*)'PROCEED WITH NEW DATA -2- ' WRITE(*,*)'PROCEED WITH OLD DATA -3- ' WRITE(*,*)'SAVE DATA FROM VIEW TO DISC -4- ' WRITE(*,*)'END -5- ' WRITE(*,'(A,$)')' ------------------ => ' READ(5,*)IAD IF (IAD.EQ.1)THEN 330 WRITE(*,*) WRITE(*,*)' NOTE: [ENTER] OR [ESC] ' WRITE(*,*)' TO ESCAPE GRAPHICS' WRITE(*,'(A,$)')' INPUT No. OF MODEL => ' READ(5,*)IVOL IF (IVOL.GT.NMOM.OR.IVOL.LT.1)GOTO 330 C !!REMOVE THIS LINE IF YOU USE OTHER COMPILER CALL GRAF(CC,CM,B,H,TEP,GRAV,NSYN,NOBS,NMOM1,NBOD,IVOL P ,AM,AL,BV,IEPER,IEPP,CA,CB,CD,CE) GOTO 320 ENDIF IF (IAD.EQ.2) GOTO 10 IF (IAD.EQ.3) GOTO 15 IF (IAD.EQ.4)THEN WRITE(*,*) WRITE(*,'(A,$)')' No. OF MODEL => P ' READ(5,*)IVOL WRITE(*,'(A,$)')' FILE NAME WITH OBSERVED PROFIL => P ' READ(5,'(A12)')INN1 WRITE(*,'(A,$)')' FILE NAME WITH SYNTHETIC PROFIL => P ' READ(5,'(A12)')INN2 CALL NAHRAJ(CC,CM,B,H,IVOL,NSYN,NOBS,NBOD,NMOM1,INN1 P ,INN2,AM,TEP,GRAV,IEPER,IEPP,BV,AL,CA,CB,CD,CE) GOTO 320 ENDIF IF (IAD.EQ.5) GOTO 400 C GOTO 10 GOTO 320 400 CLOSE(9) STOP END C----------------------------------------------------------------------- C PODPROGRAMY C----------------------------------------------------------------------- SUBROUTINE NORM(I,NMOM1,NSYN,AL,CC,H) C IT WILL NORMALIZE SYNTHETIC PROFILES IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION CC(NSYN),H(NMOM1,NSYN) IF (AL.GT.CC(NSYN))THEN AP=(1000.-H(I,NSYN))*(AL-CC(NSYN))/(10.-CC(NSYN))+H(I,NSYN) GOTO 110 ENDIF DO 100 IL=1,NSYN IF(CC(IL)-AL) 100,80,90 80 AP=H(I,IL) GOTO 110 90 AP=(H(I,IL)-H(I,IL-1))*(AL-CC(IL-1))/(CC(IL)-CC(IL-1)) P +H(I,IL-1) GOTO 110 100 CONTINUE 110 IF(ABS(AP-0.).LT.1.E-8) AP=1000. DO 120 J=1,NSYN H(I,J)=H(I,J)/AP*1000. IF (H(I,J).GT.1000.) H(I,J)=1000. 120 CONTINUE RETURN END C----------------------------------------------------------------------- FUNCTION AKONV(XX,X,F,SPACE,NSPACE,CUTOFF,NSYN2) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C X(I) - POLE X-OVYCH HODNOT, X(I+1) > X(I) C F(I) - POLE Y-OVYCH HODNOT FILTROVANEJ FUNKCIE V BODOCH X(I) C XX - X-OVA SURADNICA BODU KDE TA ZAUJIMA KONVOLUCIA C SPACE - INTEGRACNY KROK C NSPACE- POCET KROKOV V SIRKE FILTRA C CUTOFF- POLOVICA SIRKY FILTRA C NSYN2 - POCET DAT C AKONV - HODNOTA KONVOLUCIE (VYFITROVANA HODNOTA FUNKCIE) V XX C ROT - FILTRUJUCA FUNKCIA (JE VYHODNE ABY JEJ INTEGRAL=1 C => FILTROVANIM SA NEMENI INTEGRAL FILTROVANEJ FUNKCIE) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(NSYN2),F(NSYN2) SUMA=0.D0 XI1=XX-CUTOFF CALL INTEP(XI1,YI1,X,F,NSYN2,IER) F1=ROT(XI1-XX,CUTOFF,UKOEF) DO 100 I=1,NSPACE XI2=XI1+SPACE/2.D0 XI3=XI1+SPACE F2=ROT(XI2-XX,CUTOFF,UKOEF) F3=ROT(XI3-XX,CUTOFF,UKOEF) CALL INTEP(XI2,YI2,X,F,NSYN2,IER) CALL INTEP(XI3,YI3,X,F,NSYN2,IER) SUMA=(F1*YI1+4.D0*F2*YI2+F3*YI3)+SUMA YI1=YI3 XI1=XI3 F1=F3 100 CONTINUE AKONV=SPACE*SUMA/6.D0 RETURN END C----------------------------------------------------------------------- FUNCTION ROT(DL,DLL,UKOEF) C DL - ARGUMENT ROTATION FUNCTION IN [YU] C DLL - ROTATION HALFWIDTH IN [YU] C UKOEF- LIMB DARKENING COEF. C OUTPUT - ROT. PROFIL IN POINT DL IMPLICIT DOUBLE PRECISION (A-H,O-Z) IF(DL*DL.GT.DLL*DLL)THEN ROT=0.D0 RETURN ENDIF C ROT=(2.D0*(1.D0-UKOEF)*DSQRT(1.D0-DL*DL/DLL/DLL)+0.5D0*3.1415D0 C P *UKOEF*(1.D0-DL*DL/DLL/DLL))/(3.1415D0*DLL*(1.D0-UKOEF/3.D0)) ROT=(0.6366198D0*(1.D0-UKOEF)*DSQRT(1.D0-DL*DL/DLL/DLL)+0.5D0 P *UKOEF*(1.D0-DL*DL/DLL/DLL))/(DLL*(1.D0-UKOEF*0.3333333D0)) RETURN END C----------------------------------------------------------------------- SUBROUTINE STEM(IBA,IEPER,BV,UKOEF,DLL) C OUTPUT: UKOEF - LIMB DARK. COEFF. C DLL[NM] - DOPPLER WIDTH IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION G5(50),G6(50),G(50) DO 5 I=1,50 IF (IEPER.EQ.5) G(I)=G5(I) IF (IEPER.EQ.6) G(I)=G6(I) 5 CONTINUE DO 10 I=0,10 IF (ABS(G(IBA-I)-0.).LT.1.E-20)GOTO 10 IBE=IBA-I GOTO 20 10 CONTINUE 20 DO 30 I=0,10 IF (ABS(G(IBA+I)-0.).LT.1.E-20)GOTO 30 IBF=IBA+I GOTO 40 30 CONTINUE 40 IF (IBA.EQ.IBE)THEN UKOEF=G(IBA) ELSE UKOEF=FLOAT(IBA-IBE)/FLOAT(IBF-IBE)*(G(IBF)-G(IBE))+G(IBE) ENDIF WRITE(*,*) WRITE(*,*)'LIMB DARKENING COEFFICIENT',UKOEF WRITE(*,*) IF (IEPER.EQ.5)THEN DLL=BV*434.0/299792. ELSE DLL=BV*410.2/299792. ENDIF C DEFINED BY EQUATION: I=I0[1-UKOEF+UKOEF*COS(THETA)] C FROM GRYGAR ... C LIMB DARK. COEFF. -U, FOR H5------------ DATA G5/0.,0.,0.,0.,0.,0.,0.,.35,.36,.38,0.,.4,0.,.44,0.,.48 P ,0.,.52,0.,.56,0.,0.,.69,0.,0.,.75,0.,0.,0.,.78,0.,0.,0.,0. P ,0.,0.,.8,0.,0.,0.,0.,.77,0.,0.,.84,0.,0.,0.,0.,.92/ C LIMB DARK. COEFF. -U, FOR H6------------ DATA G6/0.,0.,0.,0.,0.,0.,0.,.36,.37,.39,0.,.42,0.,.46,0.,.5 P ,0.,.54,0.,.58,0.,0.,.74,0.,0.,.8,0.,0.,0.,.82,0.,0.,0.,0. P ,0.,0.,.83,0.,0.,0.,0.,.83,0.,0.,.89,0.,0.,0.,0.,.97/ RETURN END C----------------------------------------------------------------------- C SUBROUTINE I N T E P : SPLINE INTERPOLATION SCHEME BASED C ON HERMITE POLYNOMIALS C INPUT : XP - THE CHOSEN ARGUMENT VALUE C X - THE VECTOR OF INDEPENDENT VALUES C F - THE VECTOR OF FUNCTION OR DEPENDENT VALUES C N - THE NUMBER OF POINTS IN THE (X,F) VECTORS C OUTPUT : P - THE RESULTANT INTERPOLATED VALUE C IER - THE RESULTANT ERROR PARAMETER SUBROUTINE INTEP(XP,P,X,F,N,IER) IMPLICIT DOUBLE PRECISION(A-H,O-Z) REAL LP1,LP2,L1,L2 DIMENSION F(N),X(N) IER=1 IO=1 IUP=0 IF(X(2).LT.X(1)) IUP=1 N1=N-1 IF((XP.GE.X(N).AND.IUP.EQ.0).OR.(XP.LE.X(N).AND.IUP.EQ.1)) THEN P=F(N) IER=2 RETURN ELSE IF((XP.LE.X(1).AND.IUP.EQ.0).OR. * (XP.GE.X(1).AND.IUP.EQ.1)) THEN P=F(1) IER=2 RETURN END IF DO 1 I=IO,N IF(XP.LT.X(I).AND.IUP.EQ.0) GOTO 2 IF(XP.GT.X(I).AND.IUP.EQ.1) GOTO 2 1 CONTINUE P=F(N) IER=2 RETURN 2 I=I-1 IF(I.EQ.IO-1) GOTO 4 IO=I+1 LP1=1./(X(I)-X(I+1)) LP2=1./(X(I+1)-X(I)) IF(I.EQ.1) FP1=(F(2)-F(1))/(X(2)-X(1)) IF(I.EQ.1) GOTO 3 FP1=(F(I+1)-F(I-1))/(X(I+1)-X(I-1)) 3 IF(I.GE.N1) FP2=(F(N)-F(N-1))/(X(N)-X(N-1)) IF(I.GE.N1) GOTO 4 FP2=(F(I+2)-F(I))/(X(I+2)-X(I)) 4 XPI1=XP-X(I+1) XPI=XP-X(I) L1=XPI1*LP1 L2=XPI*LP2 P=F(I)*(1.-2.*LP1*XPI)*L1*L1+F(I+1)*(1.-2.*LP2*XPI1) * *L2*L2+FP2*XPI1*L2*L2+FP1*XPI*L1*L1 RETURN END C----------------------------------------------------------------------- C !!REMOVE THIS ROUTINE IF YOU USE OTHER COMPILER SUBROUTINE GRAF(CC,CM,B,H,TEP,GRAV,NSYN,NOBS,NMOM1,NBOD P ,IVOL,AM,AL,BV,IEPER,IEPP,CA,CB,CD,CE) C GRAPHICS PARAMETER (NOBS1=10000,NSYN1=26) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION CC(NSYN),CM(NOBS),B(NOBS),H(NMOM1,NSYN),TEP(NMOM1) P ,GRAV(NMOM1),IC(NSYN1),ICM(NOBS1),IB(NOBS1),IH(NSYN1) SCAL1=640./6. SCAL2=480./1200. IP1=40 IP2=480-30 DO 10 J=1,NSYN IC(J)=INT(CC(J)*SCAL1)+IP1 IH(J)=-INT(H(IVOL,J)*SCAL2)+IP2 10 CONTINUE DO 20 J=1,NBOD ICM(J)=INT(DABS(CM(J))*SCAL1)+IP1 IB(J)=-INT(B(J)*SCAL2)+IP2 20 CONTINUE CALL VGA@ C X-AXIS CALL DRAW_LINE@(IP1,IP2,640,IP2,2) DO 30 J=1,5 CALL DRAW_LINE@(INT(J*SCAL1)+IP1,IP2,INT(J*SCAL1)+IP1,IP2-5,2) 30 CONTINUE C Y-AXIS CALL DRAW_LINE@(IP1,IP2,IP1,0,2) DO 40 J=1,5 CALL DRAW_LINE@(IP1,-INT(200*J*SCAL2)+IP2,IP1+3 P ,-INT(200*J*SCAL2)+IP2,2) 40 CONTINUE C OBSERVED (CORRECTED TO SCATTERED LIGHT) PROFIL DO 50 J=1,NBOD CALL CROSS(ICM(J),IB(J),2,5) 50 CONTINUE C SYNTHETIC (CORRECTED ...) PROFIL CALL POLYLINE@(IC,IH,NSYN,1) C DESCRIPTION OF THE GRAPH CALL SET_CURSOR_POS@(47,11) WRITE(*,60)TEP(IVOL) 60 FORMAT('T=',F7.0,' K') CALL SET_CURSOR_POS@(47,13) WRITE(*,70)GRAV(IVOL) 70 FORMAT('LOG G=',F5.2,' [CGS]') CALL SET_CURSOR_POS@(47,15) WRITE(*,80)IEPER 80 FORMAT('LINE: H',I1) CALL SET_CURSOR_POS@(47,17) WRITE(*,90)10.0**FLOAT(-IEPP) 90 FORMAT('MET.ABUND./SOLAR=',F4.2) BVDEL=0. IF (TEP(IVOL).GE.CA.AND.TEP(IVOL).LE.CB.AND.GRAV(IVOL).GE.CD P .AND.GRAV(IVOL).LE.CE)BVDEL=BV CALL SET_CURSOR_POS@(47,19) WRITE(*,100)BVDEL 100 FORMAT('ROT.VEL.=',F4.0,' KM/S') CALL SET_CURSOR_POS@(47,21) WRITE(*,110)AM 110 FORMAT('SCAT.LIGHT:',F4.1,' %') CALL SET_CURSOR_POS@(47,23) WRITE(*,120)AL 120 FORMAT('NORM.:',F4.1,' NM') 130 CALL GET_KEY@(IVOL) IF (IVOL.NE.13.AND.IVOL.NE.27)GOTO 130 CALL TEXT_MODE@ RETURN END C----------------------------------------------------------------------- C !!REMOVE THIS ROUTINE IF YOU USE OTHER COMPILER SUBROUTINE CROSS(IX,IY,ILENG,ICOLL) IMPLICIT DOUBLE PRECISION(A-H,O-Z) CALL DRAW_LINE@(IX-ILENG,IY,IX+ILENG,IY,ICOLL) CALL DRAW_LINE@(IX,IY-ILENG,IX,IY+ILENG,ICOLL) RETURN END C----------------------------------------------------------------------- SUBROUTINE NAHRAJ(CC,CM,B,H,IVOL,NSYN,NOBS,NBOD,NMOM1,INN1,INN2 P ,AM,TEP,GRAV,IEPER,IEPP,BV,AL,CA,CB,CD,CE) IMPLICIT DOUBLE PRECISION(A-H,O-Z) CHARACTER*12 INN1,INN2 DIMENSION CC(NSYN),CM(NOBS),B(NOBS),H(NMOM1,NSYN), P TEP(NMOM1),GRAV(NMOM1) OPEN(6,FILE=INN1,STATUS='UNKNOWN') OPEN(7,FILE=INN2,STATUS='UNKNOWN') WRITE(6,110)AM 110 FORMAT('SCAT.LIGHT:',F4.1,' %') DO 10 I=1,NBOD WRITE(6,30)CM(I),B(I)/1000. 10 CONTINUE WRITE(7,60)TEP(IVOL) 60 FORMAT('T=',F7.0,' K') WRITE(7,70)GRAV(IVOL) 70 FORMAT('LOG G=',F5.2,' [CGS]') WRITE(7,80)IEPER 80 FORMAT('LINE: H',I1) WRITE(7,90)10.0**FLOAT(-IEPP) 90 FORMAT('MET.ABUND./SOLAR=',F4.2) BVDEL=0. IF (TEP(IVOL).GE.CA.AND.TEP(IVOL).LE.CB.AND.GRAV(IVOL).GE.CD P .AND.GRAV(IVOL).LE.CE)BVDEL=BV WRITE(7,100)BVDEL 100 FORMAT('ROT.VEL.=',F5.0,' KM/S') WRITE(7,120)AL 120 FORMAT('NORM.:',F4.1,' NM') DO 20 I=1,NSYN WRITE(7,40)CC(I),H(IVOL,I)/1000. 20 CONTINUE 30 FORMAT(F6.3,2F7.3) 40 FORMAT(F6.3,F7.3) CLOSE(6) CLOSE(7) RETURN END C----------------------------------------------------------------------- C----------------------------------------------------------------------- SUBROUTINE ZORAD(D,OMM,LL,NMOM,IVOL,NMOM1) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION D(NMOM1),OMM(NMOM1),LL(NMOM1) DO 260 I=1,IVOL OMM(I)=1.E20 LL(I)=0 260 CONTINUE DO 300 J=1,IVOL DO 290 I=1,NMOM DO 280 IL=1,IVOL IF (I.EQ.LL(IL).OR.ABS(D(I)-0.).LT.1.E-20) GOTO 290 280 CONTINUE IF (D(I).LT.OMM(J)) THEN OMM(J)=D(I) LL(J)=I ENDIF 290 CONTINUE 300 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE SYNPRO(TEP,GRAV,H,NMOM1,NMOM,NSYN,INTRP,NMOD, P IEPER,IEPP) C THIS READS SYNTHETIC PROFILES C NMOM1 - NUMBER OF MODELS(+20 INTERP.MODELS), C NSYN - NUMBER OF TABULATED POINTS IN BALM. PROFIL C NOBS - NUMBER OF OBSERVED POINTS IN BALM. PROFIL C NEZ2 - DIMENZATION PARAMETER IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION TEP(NMOM1),GRAV(NMOM1),H(NMOM1,NSYN) IF (IEPER.EQ.5)THEN IF (IEPP.EQ.0)OPEN(7,FILE='H5.DAT',STATUS='OLD') IF (IEPP.EQ.1)OPEN(7,FILE='H5-1.DAT',STATUS='OLD') IF (IEPP.EQ.2)OPEN(7,FILE='H5-2.DAT',STATUS='OLD') ELSE IF (IEPP.EQ.0)OPEN(7,FILE='H6.DAT',STATUS='OLD') IF (IEPP.EQ.1)OPEN(7,FILE='H6-1.DAT',STATUS='OLD') IF (IEPP.EQ.2)OPEN(7,FILE='H6-2.DAT',STATUS='OLD') ENDIF DO 70 I=1,NMOM1 READ(7,*,END=75)TEP(I),GRAV(I),(H(I,J),J=1,NSYN) 70 CONTINUE 75 NMOD=I-1 NMOM=NMOD+INTRP CLOSE(7) RETURN END C----------------------------------------------------------------------- C SUBROUTINE I N T E P 2 : SPLINE INTERPOLATION SCHEME BASED C ON HERMITE POLYNOMIALS C INPUT : X1A - GRID OF X'S (LENGHT M) C X2A - GRID OF Y'S (LENGHT N) C YA - 2D FUNCTION AT GRID POINTS C X1,X2-REGUESTED POINT C OUTPUT : Y - THE RESULTANT INTERPOLATED VALUE AT [X1,X2] C IER - THE RESULTANT ERROR PARAMETER SUBROUTINE INTEP2(X1A,X2A,YA,M,N,X1,X2,Y) PARAMETER (MMAX=100,NMAX=20) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION X1A(M),X2A(N),YA(M,N),YNTMP(NMAX),YMTMP(MMAX) DO 12 J=1,M DO 11 K=1,N YNTMP(K)=YA(J,K) 11 CONTINUE CALL INTEP(X2,YMTMP(J),X2A,YNTMP,N,IER) 12 CONTINUE CALL INTEP(X1,Y,X1A,YMTMP,M,IER) RETURN END