c-------------------------------------------------------------------- c Find Bar length from *CYL* table c-------------------------------------------------------------------- PARAMETER (Nmax = 5000) ! length of table DIMENSION Rad(Nmax),Surf(Nmax),SurfMajor(Nmax),Surfminor(Nmax), & Vrot(Nmax),sig_rad(Nmax),sig_z(Nmax),sig_phi(Nmax) REAL INPUT Character FileASCII*120,Analysis*3 Character File1*100,File2*100 Logical iEXIST C................................................................... C Read data and open files write(*,*)' enter the snapshot. Example: 0.6235' read(*,*) amoment c Analysis ='GDG' C Analysis ='PKD' Analysis ='AAA' c Analysis ='' c write(FileASCII,'("RESULTS/Result",".CYL.",f6.4,".dat")') c & amoment write(FileASCII,'("RESULTS/Result",a,".CYL.",f6.4,".dat")') & Analysis,amoment write (*,'(a)') FileASCII OPEN(1,FILE=FileASCII,STATUS='UNKNOWN') write(FileASCII,'("Profile.",f6.4,".dat")') & amoment OPEN(20,FILE=FileASCII,STATUS='UNKNOWN') LinesSkip = 9 CALL ReadTable(LinesSkip,aexpn,time,istep,Ndisk, & Rad,Surf,SurfMajor,Surfminor, & Vrot,sig_rad,sig_z,sig_phi, & Nlast,Analysis) c Print profile write (20,'(T3,a,T12,a,T23,a,T30,a)') & 'Aexpn','T(Gyrs)','istep','Ndisk' write (20,'(f8.4,T10,g12.4,T23,i5,T28,i8)') & aexpn,time,istep,Ndisk write (20,'(T5,a,T15,a,T25,a,T40,a)') & 'R(kpc)','Rminor','AxRatio','Surfdens' dx = 0.20 ! first approximation: large step r = 0. ! initial radius aLimit = 1.4 ! Bar should have larger axial ratio than Alimit ratio =0. do while (r.lt.12.) r = r+dx S = Value(r,SurfMajor,Rad,Nlast) rr= GetRad(S,SurfMinor,Rad,Nlast) ratio = r/rr write (20,'(f8.4,T15,f8.4,T25,f8.4,T40,g12.4)')r,rr,ratio,S enddo c Find Bar radius dx = 0.10 ! first approximation: large step r = 12.+dx ! initial radius ratio =0. do while (r.gt.0.3.and.ratio.lt.aLimit) r = r-dx S = Value(r,SurfMajor,Rad,Nlast) rr= GetRad(S,SurfMinor,Rad,Nlast) ratio = r/rr enddo r = r+1.1*dx ! improve accuracy of radius dx= dx/50. ratio =0. do while (r.gt.0.3.and.ratio.lt.aLimit) r = r-dx S = Value(r,SurfMajor,Rad,Nlast) rr= GetRad(S,SurfMinor,Rad,Nlast) ratio = r/rr enddo write(*,'(a,2f8.3,a,g12.3,a,f8.3)')' Axes=',r,rr, & ' Surface Dens=',S,' AxRatio=',r/rr c Open output file INQUIRE(file='barlength.dat',EXIST=iEXist) ! If not exists, write header write (*,*) 'Exist=',iEXist OPEN(10,FILE='barlength.dat',STATUS='UNKNOWN',access='append') If(.not.iEXIST) & write (10,'(T3,a,T12,a,T23,a,T30,a,T38,a,T50,a,T58,a,T67,3a)') & 'Aexpn','T(Gyrs)','istep','Ndisk',' MajorAxis', & 'AxRatio', 'SurfDens', & ' 1kpc: Vrot sig_z kms sig_r', & ' 5kpc: Vrot sig_z sig_r' 30 write (10,'(f8.4,T10,g12.4,T23,i5,T28,i8,2x,g12.4,f7.2,g12.4, & 2(2x,3f8.2))') & aexpn,time,istep,Ndisk,r,ratio,S, & Value(1.,Vrot,Rad,Nlast), Value(1.,Sig_z,Rad,Nlast), & Value(1.,sig_rad,Rad,Nlast), & Value(5.,Vrot,Rad,Nlast), Value(5.,Sig_z,Rad,Nlast), & Value(5.,sig_rad,Rad,Nlast) stop end c-------------------------------------------------------------------- c Read table: file 1, skip lines in the header SUBROUTINE ReadTable(LinesSkip,aexpn,time,istep,Ndisk, & Rad,Surf,SurfMajor,Surfminor, & Vrot,sig_rad,sig_z,sig_phi, & Nlast,Analysis) c-------------------------------------------------------------------- DIMENSION Rad(*),Surf(*),SurfMajor(*),Surfminor(*), & Vrot(*),sig_rad(*),sig_z(*),sig_phi(*) Character Text*100,Analysis*3 IF(Analysis.eq.'GDG'.or.Analysis.eq.'PKD')Then Do i=1,1 Read(1,'(a)')Text write (*,'(a)') Text EndDo read(1,*) time,ndisk write(*,'(g11.3,i9)') time,ndisk istep =0 aexpn =time Else ! ART Do i=1,2 Read(1,'(a)')Text write (*,'(a)') Text EndDo read(1,*) aexpn,da,time,istep,ndisk write(*,'(3g11.3,i5,i9)') aexpn,da,time,istep,ndisk EndIf Do i=4,LinesSkip Read(1,'(a)')Text write (*,'(a)') Text EndDo i =0 1 i=i+1 read(1,*,end=100,err=100)rout,Rad(i),Vrot(i), & sig_phi(i),sig_z(i),sig_rad(i), & Surf(i),d0,Vrot_bar,sig_rotbar,SurfMajor(i), & Vrot_off,sig_rotoff,SurfMinor(i),npart goto 1 100 Nlast=i-1 write (*,*) ' Read Nlines from table =',Nlast return end c-------------------------------------------------------------------- c get interpolated value at radius r c Table -- table of function c Rad -- table of radius FUNCTION Value(r,Table,Rad,Nlast) c-------------------------------------------------------------------- DIMENSION Rad(*),Table(*) If(r.ge.Rad(Nlast))Then Value = Table(Nlast) return EndIf If(r.le.Rad(1))Then Value = Table(1) return EndIf do i=2,Nlast if(r.le.Rad(i))Then Value =Table(i)+(Table(i-1)-Table(i))/ & (Rad(i)-Rad(i-1))*(Rad(i)-r) goto 10 EndIf EndDo write (*,*) ' error in function Value:',r,Nlast stop 10 return end c-------------------------------------------------------------------- c get interpolated value of radius r for S c Table -- table of function c Rad -- table of radius FUNCTION GetRad(S,Table,Rad,Nlast) c-------------------------------------------------------------------- DIMENSION Rad(*),Table(*) If(S.le.Table(Nlast))Then GetRad = Rad(Nlast) return EndIf do i=Nlast,2,-1 if(S.ge.Table(i).and.S.lt.Table(i-1))Then GetRad =Rad(i)-(Rad(i)-Rad(i-1))*(S-Table(i))/ & (Table(i-1)-Table(i)) goto 10 EndIf EndDo c write (*,*) ' error in function GetRad:',S,Nlast GetRad =Rad(1) 10 return end