program pktnor4 C ============================================== C ver. 3-Mar-1996 S. Zola C popr. 13-April-2003 C ver. 9-May-2003 also errors in flux added C ver.15-Jan-2005 avoids flux errors =0.000 C ============================================== parameter (mdim=5000) implicit double precision (a-h,o-z) double precision magb,magv dimension phasB(mdim),magB(mdim),fazaB(mdim) 1,phasV(mdim),magV(mdim),fazaV(mdim) character *10 fil1,fil2 common /elem/xm0,p0 common /param/x1,x2,x3,x4,s1,s2,fac common /f/factb,factv common /kolor/alfab,alfav,ss1,ss2 common /range/ r1,r2,r3,r4 common /odr/ px1,px2 common /sig/czy,stala C C data file C Please give only up to 5 digits of JDhel or else this program will not C work i.e. 50282.272727 write(*,*) ' INPUT FILE: ' read(*,'(a)') fil1 write(*,*) ' OUTPUT FILE: ' read(*,'(a)') fil2 open(1,file=fil1) open(2,file=fil2) write(*,*) ' Threshold level for removing bad points?: Nm,Nd ' C ============================================================================================================================= C These parameters should be set by trial, when seen how many C single points is being removed from anormal point. C Nm - used while removing bad points near the minima C Nd - used while removing bad points outside the minima C They could be the same, though, if increase/decrease is fast, C I suggest they are to be different to account for fast brightness C change. Please note, removing bad pints is done only once! C ============================================================================================================================= read(*,*) px1,px2 C ============================================================================================================================= C czy -> factor1 = +-1. (To change sign of magnitudes if needed) C stala -> factor2 (To subtract a constant if delta_magnitudes are big) C ============================================================================================================================= read(1,*)czy,stala write(2,*)' factor =',czy,' subtracted: ', stala C C ============================================================================================================================= C x1-x4 parameters describing phase ranges, s1,s2 - steps in the ranges, where: C between <0-x1> i normal points are calculated using step "s1" (primary minimum) C between i normal points are calculated using step "s2" C between using step "s1*fac" (secondary minimum) C ============================================================================================================================= read(1,*)x1,x2,x3,x4,s1,s2,fac if(x1.lt.-1.)stop C ============================================================================================================================= C coefficients for standard reduction C this is not done right, it is just a "zero-order" aproximation C if the coefficients are set to 0 nothing is done (recommended) C ============================================================================================================================= read(1,*)alfab,alfav read(1,*)r1,r2,r3,r4 C ------------------------------------------------------------------- C Epoch (primary minimum) Period C ------------------------------------------------------------------- read(1,*)xm0,p0 C ------------------------------------------------------------------- C === factors to normalize LC in flux units for required phase C The code must be run twice, set the coefficients to 1.0 for the initial C run and take the coefficients from the output file for the next run C ------------------------------------------------------------------- read(1,*)factb,factv write(2,*) ' ' write(2,*)' M0 P ' write(2,'(5x,f12.5,3x,f12.6)') xm0,p0 write(2,*) ' ' call czyt(phasb,magb,nobsb) call czyt(phasv,magv,nobsv) call pfaz(phasb,nobsb,fazab) call pfaz(phasv,nobsv,fazav) call sort3(fazab,phasb,magb,nobsb) call sort3(fazav,phasv,magv,nobsv) IF(ALFAb.ne.0..OR.ALFAV.NE.0.) THEN CALL systm(fazab,magb,nobsb,fazav,magv,nobsv) END IF write(2,*) ' ********************* ' do 450 i=1,nobsb pr=(10.**(-magb(i)/2.5))/factb 450 write(2,460) i,phasb(i),fazab(i),magb(i),pr write(2,*) ' ********************* ' do 451 i=1,nobsv pr=(10.**(-magv(i)/2.5))/factv 451 write(2,460) i,phasv(i),fazav(i),magv(i),pr write(2,*) ' ********************* ' call sort(fazab,magb,nobsb,1) call sort(fazav,magv,nobsv,2) stop 460 format(1x,i4,2x,f12.5,3x,f8.5,5x,f6.3,3x,f6.4) 300 format(3x,f12.4,2x,f8.4,2x,f7.3) end subroutine czyt(phas,mag,nobs) C C Reading data - give only up to 5 digits of the JDhel, otherwise C the program will not work C parameter (mdim=5000) implicit double precision (a-h,o-z) double precision mag dimension phas(mdim),mag(mdim) common /sig/czy,stala do 74 i=1,5000 read(1,*)phas(i),mag(i) mag(i)=(mag(i)*czy) + stala if(phas(i).lt.0.)stop if(phas(i)-60000.) 74,125,125 74 continue 125 nobs=i-1 write(2,30) nobs 30 format(/' liczba obserwacji =',i6/) return end subroutine pfaz(phas,nobs,faza) parameter (nbs=5000) implicit double precision (a-h,o-z) C calculating phases with given M0 and P dimension phas(nbs),faza(nbs) common /elem/ xm0,p0 do 100 i=1,nobs faza(i)=(phas(i)-xm0)/p0 t=aint(faza(i)) faza(i)=faza(i)-t if(faza(i).lt.0.) faza(i)=faza(i)+1. 100 continue return end subroutine sort(faza,mag,nobs,ik) parameter (mdim=5000,ndm=5000) implicit double precision (a-h,o-z) double precision mag dimension xfazn(ndm),xmagn(ndm),xvar(ndm),numb(ndm) dimension xfaznfl(ndm),xfln(ndm),xvarfl(ndm),numbfl(ndm) dimension mag(nobs),faza(nobs),flux(nobs) dimension fl(mdim),weight(mdim) common /f/factb,factv factor=factb if(ik.eq.2)factor=factv call norm(faza,mag,nobs,xfazn,xmagn,numb,xvar,il) il=il-1 write(2,55) do 501 i=1,il fl(i)=(10.**(-xmagn(i)/2.5))/factor C ---------------------------------------------- C Next record assumes the weight for normal points for the C DC W-D code (here the same weight = 1. but can be set C in a required way) C ---------------------------------------------- weight(i)=1. write(2,54) xfazn(i),fl(i),xvar(i),xmagn(i),numb(i),i 501 continue write(2,51) do 57 i=1,il,5 write(2,56)xfazn(i),fl(i),weight(i),xfazn(i+1),fl(i+1),weight(i+1) $,xfazn(i+2),fl(i+2),weight(i+2),xfazn(i+3),fl(i+3),weight(i+3) $,xfazn(i+4),fl(i+4),weight(i+4) 57 continue C ===================================================================== C flux errors for testing only for now C howerver, it seems to work OK C ===================================================================== do 502 i=1,nobs C flux(i)=(10.**(-mag(i)/2.5)/factor) flux(i)=10.**(-mag(i)/2.5) 502 continue call norm(faza,flux,nobs,xfaznfl,xfln,numbfl,xvarfl,iln) iln=iln-1 write(2,*) ' ' write(2,*) ' ============= FLUX errors ==================' do 503 i=1,iln if(xvarfl(i).lt.0.001) xvarfl(i)=0.001 write(2,59) xfaznfl(i),xfln(i)/factor,xvarfl(i),numbfl(i),i 503 continue write(2,*) ' ============= END FLUX errors ==============' write(2,*) ' ' C ===================================================================== 52 format(i5) 55 format(1x,'NORMAL POINTS (PHASE-FLUX-----MAGN.StDv.--NO of PTS'/) 54 format(1x,f8.4,2x,f6.4,2x,f5.3,2x,f7.3,4x,i4,I5) 51 format(1x,'N',4X) 56 format(5(f6.4,1x,f6.4,f3.1)) 59 format(1x,f8.4,2x,f6.4,2x,f5.3,4x,i4,I5) RETURN END subroutine norm(faza,fmag,nobs,xfazn,xmagn,numb,xvar,iik) parameter (ndm=5000,nbs=5000) implicit double precision (a-h,o-z) dimension faza(nbs),fmag(nbs),xfazn(ndm),xmagn(ndm) dimension xfaz(ndm),xmag(ndm),numb(ndm),xvar(ndm) common /param/x1,x2,x3,x4,s1,s2,fac x0=0. iik=1 888 continue if(x0.lt.x1.or.x0.ge.x4) then step=s1 else if(x0.gt.x2.and.x0.lt.x3)then step=s1*fac else step=s2 end if end if xk=x0+step if(xk.gt.1.05) go to 130 k=0 do 100 i=1,nobs if(faza(i).ge.x0.and.faza(i).lt.xk) go to 10 go to 100 10 k=k+1 xfaz(k)=faza(i) xmag(k)=fmag(i) 100 continue if(k.eq.0.) go to 999 call sre(xfaz,k,xf) call sre(xmag,k,xm) if(k.eq.1) xvar(iik)=10. if(k.gt.1) call vari(xmag,k,var,xm) if(k.gt.2) call odrz(xfaz,xmag,k,var,xm,xf) if(k.gt.1) xvar(iik)=var if(k.eq.1) xvar(iik)=10. xfazn(iik)=xf xmagn(iik)=xm numb(iik)=k iik=iik+1 999 x0=xk go to 888 130 return end subroutine vari(xmag,k,var,xm) C computing standard deviation parameter (ndm=5000) implicit double precision (a-h,o-z) dimension xmag(ndm),del(ndm) do 50 i=1,k 50 del(i)=xmag(i)-xm d0=0. do 70 i=1,k d1=d0+del(i)*del(i) d0=d1 70 continue var=d0/(float(k)*float(k-1)) var=sqrt(var) return end subroutine systm(xb,yb,nb,xv,yv,nv) parameter (mdim=5000) implicit double precision (a-h,o-z) C ========================================================== C aproximate reduction to the standard system C WARNING!!! IT may not work properly C alfa - coefficient C r1-r4 ranges C ========================================================== dimension xb(nb),yb(nb),xv(nv),yv(nv),color(mdim) common /param/p1,p2,p3,p4,s1,s2,fac common /kolor/alfab,alfav,sss1,sss2 common /rang/r1,r2,r3,r4 iil=1 x0=0. 999 continue if(x0.lt.r1.or.x0.ge.r4) then step=s1 if(x0.gt.r2.and.x0.lt.r3) then step=fac*s1 else step=s2 end if end if xk=x0+step if(xk.gt.1.06) return call mean(xb,yb,nb,sr1,x0,xk,k1) call mean(xv,yv,nv,sr2,x0,xk,k2) color(iil)=sr1-sr2 if(k1.eq.0.or.k2.eq.0) color(iil)=color(iil-1) do 101 i=1,nb if(xb(i).ge.x0.and.xb(i).lt.xk)yb(i)=yb(i)+color(iil)*alfab 101 continue do 102 i=1,nv if(xv(i).ge.x0.and.xv(i).lt.xk)yv(i)=yv(i)+color(iil)*alfav 102 continue 104 x0=xk iil=iil+1 go to 999 C If there is a warning in your compiler about the next line, C just ignore it return end subroutine mean(p,fm,n,s,x0,xk,k) parameter (ndm=5000) implicit double precision (a-h,o-z) dimension p(n),fm(n),xmag(ndm) k=0 do 100 i=1,n if(p(i).ge.x0.and.p(i).lt.xk) go to 10 go to 100 10 k=k+1 xmag(k)=fm(i) 100 continue if(k.eq.0) go to 350 call sre(xmag,k,s) return 350 s=30. return end subroutine sortr(x,y,n) parameter (ndm=5000) implicit double precision (a-h,o-z) C sorting dimension x(5000),y(5000) c do 10 kk=1,n do 20 i=2,n xlo=dmin1(x(i-1),x(i)) if(xlo.eq.x(i)) then tt=x(i-1) x(i-1)=xlo t=y(i-1) y(i-1)=y(i) x(i)=tt y(i)=t end if 20 continue 10 continue return end subroutine odrz(x,y,k,var,ysr,xf) C removing bad points parameter (ndm=5000) implicit double precision (a-h,o-z) dimension x(k),y(k),x1(ndm),y1(ndm) common /param/p1,p2,p3,p4,s1,s2,fac common /odr/ px1,px2 if(k.lt.2)return px=px1 if(x(1).lt.p1.or.x(1).gt.p4) px=px2 if(x(1).gt.p2.and.x(1).lt.p3) px=px2 par=px*var l=1 do 10 i=1,k p=abs(y(i)-ysr) if(p.gt.par) go to 9 x1(l)=x(i) y1(l)=y(i) l=l+1 go to 10 9 write(2,20)x(i),y(i),p,var,ysr 10 continue k1=l-1 k=k1 if(k1.eq.k) return if(k.gt.1) then call sre(x1,k1,xf) call sre(y1,k1,ysr) call vari(y1,k1,var,ysr) else xf=x1(1) ysr=y1(1) var=30. end if 20 format(1x,' pkt. wyrzucony: ',5(f8.4,1x)) return end subroutine sre(x,n,sr) parameter (ndm=5000) implicit double precision (a-h,o-z) C mean value od one dimensional matrix dimension x(ndm) a0=0. do 200 i=1,n a1=a0+x(i) a0=a1 200 continue sr=a0/float(n) return end subroutine sort3(x,y,z,n) parameter (ndm=5000) implicit double precision (a-h,o-z) dimension x(n),y(n),z(n) c do 10 kk=1,n do 20 i=2,n xlo=dmin1(x(i-1),x(i)) if(xlo.eq.x(i)) then tt=x(i-1) x(i-1)=xlo t=y(i-1) ttt=z(i-1) y(i-1)=y(i) z(i-1)=z(i) x(i)=tt y(i)=t z(i)=ttt end if 20 continue 10 continue return end