C C PROGRAM NA SKLADANIE VIACERYCH SPEKTIER (NAPR. DVOJHVIEZD), C NA CO-ADDOVANIE=SCITAVANIE X,Y DAT S MOZNYM POSUNOM C A NORMALIZACIOU ALEBO C NA JEDNODUCHU ARITMETIKU MEDZI DVOMI SPEKRAMI -X,Y SUBORMI c c authors: J. Budaj c contact: budaj@ta3.sk c published: http://www.ta3.sk/~budaj/software c special requirements: c non original routins: intep (Press.et.al. 1992, Numerical Recipes) c input: (3)-specsum.in, (2)-files c output: (4)-specsum.out c PARAMETER(NBOD=10000) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*30 IN,INN DIMENSION X(NBOD),XI(NBOD),F(NBOD),FI(NBOD),IN(100), P CONT(100),SIFT(100),NDAT(100) OPEN(3,FILE='specsum.in',STATUS='old') c WRITE(*,*)' VYBER OPERACIU: ' c WRITE(*,*)' ZAKL.SUBOR + DRUHY SUBOR [1]' c WRITE(*,*)' ZAKL.SUBOR - DRUHY SUBOR [2]' c WRITE(*,*)' ZAKL.SUBOR * DRUHY SUBOR [3]' c WRITE(*,*)' ZAKL.SUBOR / DRUHY SUBOR [4]' c WRITE(*,*)' COADDOVANIE VIACERYCH SUBOROV A NORMALIZACIA [5]' c WRITE(*,'(A,$)')' OPERACIA CISLO => ' READ(3,'(//////)') READ(3,*)NOPER WRITE(*,*)noper c WRITE(*,'(A,$)')' POCET SPEKTIER-SUBOROV => ' READ(3,*) READ(3,*)NFILE WRITE(*,*)nfile C WRITE(*,'(A,$)')' MENO ZAKL. VSTUPNEHO SUBORU => ' C READ(*,'(A30)')IN(1) C READ(3,'(A)')IN(1) C IF(NOPER.EQ.5)THEN C WRITE(*,'(A,$)')' KONTINUUM ZAKL. SPEKTRA-VAHA SUBORU => ' C READ(*,*)CONT(1) C SNORM=CONT(1) C SIFT(1)=0.D0 C ENDIF READ(3,'(//)') SNORM=0. DO 5 I=1,NFILE C WRITE(*,'(A,$)')' MENO DALSIEHO VSTUPNEHO SUBORU => ' READ(3,'(A)')IN(I) WRITE(*,'(A)')in(i) IF(NOPER.EQ.5)THEN C WRITE(*,'(A,$)')' KONTINUUM SPEKTRA-VAHA SUBORU => ' C WRITE(*,'(A,$)')' SIFT SPEKTRA V [TJ]-O KOLKO HO POSUNUT? => ' READ(3,*)CONT(I),SIFT(I) WRITE(*,*)CONT(I),SIFT(I) SNORM=SNORM+CONT(I) ELSE READ(3,*) ENDIF 5 CONTINUE c WRITE(*,*) c WRITE(*,'(A,$)')' MENO VYSTUPNEHO SUBORU => ' c READ(*,'(A30)')INN inn='specsum.out' OPEN(4,FILE=INN,STATUS='unknown') C C PRE NOPER=5 C CLIGHT=2.99792458E5 IF(NOPER.EQ.5)THEN DO 60 J=1,NFILE OPEN(2,FILE=IN(J),STATUS='OLD') I=1 10 READ(2,*,END=20)XI(I),FI(I) IF(J.EQ.1)THEN X(I)=XI(I)+SIFT(J)/CLIGHT*XI(I) F(I)=0.D0 ENDIF I=I+1 GOTO 10 20 NDAT(J)=I-1 DO 50 I=1,NDAT(1) IF(J.GT.1)THEN C VRATI POSUNUTU INTERPOLOVANU Y-HODNOTU 'POSUN' POSUNX=X(I)-SIFT(J)/CLIGHT*X(I) C CALL INTEP(X(I)-SIFT(J),POSUN,XI,FI,NDAT(J),IER) CALL INTEP(POSUNX,POSUN,XI,FI,NDAT(J),IER) ELSE POSUN=FI(I) ENDIF F(I)=CONT(J)/SNORM*POSUN+F(I) IF(J.EQ.NFILE)WRITE(4,'(2E15.7)')X(I),F(I) 50 CONTINUE CLOSE (2) 60 CONTINUE GOTO 200 ENDIF C C PRE NOPER= 1,2,3,4 C DO 160 J=1,NFILE OPEN(2,FILE=IN(J),STATUS='OLD') I=1 110 READ(2,*,END=120)XI(I),FI(I) IF(J.EQ.1)THEN X(I)=XI(I) F(I)=FI(I) ENDIF I=I+1 GOTO 110 120 NDAT(J)=I-1 IF(J.GT.1)THEN DO 150 I=1,NDAT(1) CALL INTEP(X(I),POSUN,XI,FI,NDAT(J),IER) IF(NOPER.EQ.1)F(I)=F(I)+POSUN IF(NOPER.EQ.2)F(I)=F(I)-POSUN IF(NOPER.EQ.3)F(I)=F(I)*POSUN IF(NOPER.EQ.4)F(I)=F(I)/POSUN WRITE(4,'(2E15.7)')X(I),F(I) 150 CONTINUE ENDIF CLOSE (2) 160 CONTINUE 200 CLOSE (4) STOP 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*8 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