program phical integer dim parameter (dim=200) real radius,window real vv(dim),hh(dim),vh(dim),hv(dim),phi(dim) real vva(dim),hha(dim),vha(dim),hva(dim) real vvd,hhd,vhd,hvd real sigma(200),var(4,4) real avv(4),ahh(4),avh(4),ahv(4) real chi,deg,basis(4) real mvv,mhh,mvh,mhv real vvt(dim),hht(dim),vht(dim),hvt(dim) integer ntrace,mark,delta,lb,ub,number,ma,list(4) integer ncvm,fit character filename*13,label*11,freq*2 external fun c cccccccccccccccccccccc c open file #1 c cccccccccccccccccccccc print *, 'Enter file name' read(*,'(a13)') filename c print *, 'Enter arch radius' c read *, radius open(unit=20,file=filename) read(20,*) ntrace do 10 i=1,ntrace read(20,*) vv(i),hh(i),vh(i),hv(i) vv(i)=10.0**(vv(i)/10.) hh(i)=10.0**(hh(i)/10.) vh(i)=10.0**(vh(i)/10.) hv(i)=10.0**(hv(i)/10.) phidelta=120.0/(ntrace-1) phi(i)=abs(60.0-(i-1)*phidelta) c phi(i)=abs(45.0-(i-1)*phidelta) 10 continue close(20) mark=ntrace goto 40 c cccccccccccccccccccccc c open file #2 c cccccccccccccccccccccc print *, 'Enter file name' read *, filename if(filename .eq. 'stop') then goto 40 else continue endif print *, 'Enter arch radius' read *, radius open(unit=20,file=filename) read(20,*) ntrace do 20 i=mark+1,mark+ntrace,1 read(20,*) vv(i),hh(i),vh(i),hv(i) vv(i)=10.0**(vv(i)/10.) hh(i)=10.0**(vv(i)/10.) vh(i)=10.0**(vv(i)/10.) hv(i)=10.0**(vv(i)/10.) c phidelta= c phi(i)= 20 continue close(20) mark=ntrace+mark c cccccccccccccccccccccc c open file #3 c cccccccccccccccccccccc print *, 'Enter file name' read *, filename if(filename .eq. 'stop') then goto 40 else continue endif print *, 'Enter arch radius' read *, radius open(unit=20,file=filename) read(20,*) ntrace do 30 i=mark+1,mark+ntrace,1 read(20,*) vv(i),hh(i),vh(i),hv(i) vv(i)=10.0**(vv(i)/10.) hh(i)=10.0**(vv(i)/10.) vh(i)=10.0**(vv(i)/10.) hv(i)=10.0**(vv(i)/10.) c phi(i)= 30 continue close(20) mark=ntrace+mark 40 continue c cccccccccccccccccccccc c open save files c cccccccccccccccccccccc label=filename(3:13) freq=filename(1:2) open(unit=21,file=freq//'d'//label) open(unit=22,file=freq//'a'//label) open(unit=23,file=freq//'t'//label) open(unit=24,file=freq//'c'//label) c cccccccccccccccccccccc c sort data in phi c cccccccccccccccccccccc call piksrt(mark,phi,vv,hh,vh,hv) c cccccccccccccccccccccc c smooth data using c moving average c cccccccccccccccccccccc delta=10 window=3.5 mvv=0.0 mhh=0.0 mhv=0.0 mvh=0.0 data vva /200 * 0.0/, hha /200 * 0.0/, * hva /200 * 0.0/, vha /200 * 0.0/ do 50 i= 1,mark if (i-delta .lt. 1) then lb=1 else lb=i-delta endif if (i+delta .gt. mark) then ub=mark else ub=i+delta endif number=0 do 60 j=lb,ub if(abs(phi(i)-phi(j)) .lt. window)then vva(i)=vva(i)+vv(j) hha(i)=hha(i)+hh(j) vha(i)=vha(i)+vh(j) hva(i)=hva(i)+hv(j) number=number+1 else continue endif 60 continue vvd=10.0*alog10(vv(i)) hhd=10.0*alog10(hh(i)) vhd=10.0*alog10(vh(i)) hvd=10.0*alog10(hv(i)) vva(i)=10.0*alog10(vva(i)/number) hha(i)=10.0*alog10(hha(i)/number) vha(i)=10.0*alog10(vha(i)/number) hva(i)=10.0*alog10(hva(i)/number) c cccccccccccccccccccccc c save smoothed and c unsmoothed data c cccccccccccccccccccccc write(21,*)phi(i),vvd,hhd,vhd,hvd write(22,*)phi(i),vva(i),hha(i),vha(i),hva(i) c cccccccccccccccccccccc c determine averages c c cccccccccccccccccccccc mvvd=mvvd+vvd/mark mhhd=mhhd+hhd/mark mvhd=mvhd+vhd/mark mhvd=mhvd+hvd/mark 50 continue c cccccccccccccccccccccc c determine least squared c model for data c cccccccccccccccccccccc data avv /4 * 0.0/ data ahh /4 * 0.0/ data avh /4 * 0.0/ data ahv /4 * 0.0/ data sigma /200 * 1.0/ ma=4 list(1)=1 list(2)=2 list(3)=4 ncvm=3 fit=3 call lfit(phi,vva,sigma,mark,avv,ma,list,fit,var,ncvm,chi,fun) write(24,*) avv(1),avv(2),avv(3),avv(4),mvv data sigma /200 * 1.0/ ma=4 list(1)=1 list(2)=2 list(3)=4 ncvm=3 fit=3 call lfit(phi,hha,sigma,mark,ahh,ma,list,fit,var,ncvm,chi,fun) write(24,*) ahh(1),ahh(2),ahh(3),ahh(4),mhh data sigma /200 * 1.0/ ma=4 list(1)=1 list(2)=2 list(3)=3 ncvm=3 fit=3 call lfit(phi,vha,sigma,mark,avh,ma,list,fit,var,ncvm,chi,fun) write(24,*) avh(1),avh(2),avh(3),avh(4),mvh data sigma /200 * 1.0/ ma=4 list(1)=1 list(2)=2 list(3)=3 ncvm=3 fit=3 call lfit(phi,hva,sigma,mark,ahv,ma,list,fit,var,ncvm,chi,fun) write(24,*) ahv(1),ahv(2),ahv(3),ahv(4),mhv c cccccccccccccccccccccc c generate and store c model data c cccccccccccccccccccccc do 70 i=1,31 deg=(i-1)*2.0 call fun(deg,basis,ma) vvt(i)=avv(1)*basis(1)+avv(2)*basis(2) * +avv(3)*basis(3)+avv(4)*basis(4) hht(i)=ahh(1)*basis(1)+ahh(2)*basis(2) * +ahh(3)*basis(3)+ahh(4)*basis(4) vht(i)=avh(1)*basis(1)+avh(2)*basis(2) * +avh(3)*basis(3)+avh(4)*basis(4) hvt(i)=ahv(1)*basis(1)+ahv(2)*basis(2) * +ahv(3)*basis(3)+ahv(4)*basis(4) write(23,*) deg,vvt(i),hht(i),vht(i),hvt(i) 70 continue close(21) close(22) close(23) stop end subroutine fun(xx,afun,num) real xx,afun(4),pi integer num pi=4*atan(1.0) num=4 afun(1)=1.0 afun(2)=cos(xx*pi/180.) afun(3)=cos(8*xx*pi/180.) afun(4)=sin(xx*pi/180.) return end