inv_taup.ps

来自「著名的seismiclab的代码 是地震学研究人员必备的工具」· PS 代码 · 共 365 行

%!PS-Adobe-3.0
%%Creator: A2ps version 4.3
%%CreationDate: Thu Oct 18 21:50:12 2001
%%Pages: (atend)
%%DocumentFonts: Courier Courier-Bold Helvetica Helvetica-Bold
%%EndComments
% Copyright (c) 1993, 1994, Miguel Santana, [e-mail]

/$a2psdict 100 dict def
$a2psdict begin

% General macros.
/xdef {exch def} bind def
/getfont {exch findfont exch scalefont} bind def

% Create Courier backspace font
/backspacefont {
    /Courier findfont dup length dict begin
	{ % forall
	    1 index /FID eq { pop pop } { def } ifelse
	} forall
	currentdict /UniqueID known { % if
	    /UniqueID UniqueID 16#800000 xor def
	} if
	CharStrings length 1 add dict begin
	    CharStrings { def } forall
	    /backspace { -600 0 0 0 0 0 setcachedevice } bind def
	    currentdict
	end
	/CharStrings exch def
	/Encoding Encoding 256 array copy def
	Encoding 8 /backspace put
	currentdict
    end
    definefont pop
} bind def

% FUNCTIONS

% Function filename: Initialize file printing.
/fn
{ /filenm xdef
  /filenmwidth filenm stringwidth pop def
  /filenmfont
       filenmwidth fns gt
       {
	       filenmfontname
	       fnfs fns mul filenmwidth div
	     getfont
       }
       { sfnf }
     ifelse
  def
} bind def

% Function header: prints page header. no page
% is passed as argument.
/hp
  { x sd get  y sd get hs sub 1 add  moveto
    df setfont
    gsave
      x sd get y sd get moveto
      0 hs 2 div neg rmoveto 
      hs setlinewidth
      0.95 setgray
      pw 0 rlineto stroke
    grestore
    gsave
      dfs hm rmoveto
      d show				% date/hour
    grestore
    gsave
      pnum cvs pop				% page pop up
        pw (Page 999) stringwidth pop sub
        hm
	  rmoveto
      (Page ) show pnum show		% page number
    grestore
    empty pnum copy pop
    gsave
      filenmfont setfont
         fns filenm stringwidth pop sub 2 div dw add
          bm 2 mul 
        add 
        hm
      rmoveto
        filenm show			% file name
      grestore
    } bind def

% Function border: prints border page
/border 
{ x sd get y sd get moveto
  gsave				% print four sides
    0.7 setlinewidth		% of the square
    pw 0 rlineto
    0 ph neg rlineto
    pw neg 0 rlineto
    closepath stroke
  grestore
} bind def

% Function hborder: completes border of the header.
/hborder 
{ gsave
	0.7 setlinewidth
	0 hs neg rmoveto
	pw 0 rlineto
	stroke
  grestore
} bind def

% Function sheetnumber: prints the sheet number.
/sn
    { snx sny moveto
      df setfont
      pnum cvs
	  dup stringwidth pop (0) stringwidth pop sub neg 0 rmoveto show
      empty pnum copy pop
    } bind def

% Function loginprint: prints the login id of the requestor.
/lgp
    { lx ly moveto
      df setfont
      dup stringwidth pop neg 0 rmoveto show
    } bind def

% Function currentdate: prints the current date.
/cd
    { dx dy moveto
      df setfont
      (Printed: ) show
      td show
    } bind def

% Function filename_footer: prints the file name at bottom of page.
/fnf
    { fnx fny moveto
      df setfont
      filenm center show
    } bind def

% Function center: centers text.
/center
    { dup stringwidth pop
      2 div neg 0 rmoveto
    } bind def

% Function s: print a source line
/s  { show
      /y0 y0 bfs sub def
      x0 y0 moveto
    } bind def

% Functions b and st: change to bold or standard font
/b  { show
      bdf setfont
    } bind def
/st { show
      bf setfont
    } bind def

% Strings used to make easy printing numbers
/pnum 12 string def
/empty 12 string def

% Global initializations

/CourierBack backspacefont
/filenmfontname /Helvetica-Bold def
/inch {72 mul} bind def

%
% Meaning of some variables and functions (coded names)
%
%	twp:		twinpages?
%	sd:		sheet side
%	l:		line counter
%	c:		column counter
%	d:		date
%	td:		current date (for today)
%	lg:		login name
%	fn:		filename printing function
%	sn:		sheetnumber printing function
%	cd:		current date printing function
%	fnf:		filename footer printing function
%	lgp:		login printing function
%	hp:		header printing function
%	y:		y coordinate for the logical page
%	x:		x coordinate for the logical page
%	sny:		y coordinate for the sheet number
%	snx:		x coordinate for the sheet number
%	dy:		y coordinate for the date
%	dx:		x coordinate for the date
%	ly:		y coordinate for the login
%	lx:		x coordinate for the login
%	scx:		x coordinate for the sheet center
%	fny:		y coordinate for the filename (footer)
%	fnx:		x coordinate for the filename (footer)
%	fnfs:		filename font size
%	bfs:		body font size
%	dfs:		date font size
%	bfs:		body font size
%	df:		date font
%	bf:		body font
%	bdf:		bold font
%	sfnf:		standard filename font
%	dw:		date width
%	pw:		page width
%	sw:		sheet width
%	ph:		page height
%	sh:		sheet height
%	hm:		header margin
%	tm:		top margin
%	bm:		body margin
%	rm:		right margin
%	lm:		left margin
%	hs:		header size
%	fns:		filename size

% Initialize page description variables.
/x0 0 def
/y0 0 def
/sh 11 inch def
/sw 8.5 inch def
/margin 1 inch def
/rm margin 3 div def
/lm margin 2 mul 3 div def
/d () def
/td (Oct 18 2001 21:50) def
/lg (Printed from bayes.phys.ualberta.ca on) def
%%EndProlog

/docsave save def
%%Page: 1 1
/pagesave save def
/twp true def
/fnfs 11 def
/dfs fnfs 0.8 mul def
/df /Helvetica dfs getfont def
/dw df setfont td stringwidth pop def
/sfnf filenmfontname fnfs getfont def
/hm fnfs 0.25 mul def
/hs 0.22 inch def
/bfs 6.8 def
/bdf /Courier-Bold bfs getfont def
/bm bfs 0.7 mul def
/bf /CourierBack bfs getfont def
/l 76 def
/c 87 def
/pw
   bf setfont (0) stringwidth pop c mul bm dup add add
   def
/ph
   bfs l mul bm dup add add hs add
   def
/fns
      pw
      fnfs 4 mul dw add (Page 999) stringwidth pop add
    sub
  def
/tm margin twp {3} {2} ifelse div def
/sd 0 def
/y [ rm ph add bm add
          dup ] def
/sny dfs dfs add def
/snx sh tm dfs add sub def
/dy sny def
/dx tm dfs add def
/x [ tm			% left page
          dup 2 mul pw add	% right page
        ] def
/scx sh 2 div def
/fny dy def
/fnx scx def
/ly fnfs 2 div y sd get add def
/lx snx def
/d (Sep 29 2001 00:17) def
( inv_taup.m ) fn
sw 0 translate
90 rotate
1 hp
border
hborder
/x0 x 0 get bm add def
/y0 y 0 get bm bfs add hs add sub def
x0 y0 moveto
bf setfont
( function [m] = inv_taup\(d,dt,h,q,N,flow,fhigh,mu\);) s
( %INV_TAUP: Inverse Radon transform) s
( %) s
( %  Given the seismic data this subroutine is used) s
( %  to compute the  Radon panel. By Radon Transform) s
( %  we mean the Linear or Parabolic Radon Transform.) s
( %) s
( %  [m] = inv_taup\(d,dt,h,q,N,flow,fhigh,mu\)) s
( % ) s
( %  IN   d: seismic traces  \(d\(nt,nh\) ) s
( %       dt: sampling in sec) s
( %       h\(nh\) offset or position of traces in mts) s
( %       q\(nq\) ray parameters to retrieve or curvature) s
( %             of the parabola if N=2) s
( %       N:1 Linear tau-p  ) s
( %        :2 Parabolic tau-p) s
( %       flow:  freq.  where the inversion starts in HZ \(> 0Hz\)) s
( %       fhigh: freq. where the inversion ends in HZ   \(< Nyquist\) ) s
( %       mu: regularization parameter ) s
( %) s
( %  OUT  m: the linear or parabolic tau-p panel [m\(nt,nq\)]) s
( %) s
( %  SeismicLab) s
( %  Version 1) s
( %) s
( %  written by M.D.Sacchi, last modified December 10,  1998.) s
( %  [e-mail]) s
( %) s
( %  Copyright \(C\) 1998 Signal Analysis and Imaging Group) s
( %                     Department of Physics) s
( %                     The University of Alberta) s
( %) s
( ) s
( nt= max\(size\(d\)\);) s
( nq = max\(size\(q\)\);) s
( nh = max\(size\(h\)\);) s
( ) s
( D = fft\(d,[],1\);) s
( M = zeros\(nt,nq\);) s
( i = sqrt\(-1\);) s
( ) s
( ilow  = floor\(flow*dt*nt\)+1; if ilow<1; ilow=1;end;) s
( ihigh = floor\(fhigh*dt*nt\)+1;) s
( if ihigh>floor\(nt/2\)+1; ihigh=floor\(nt/2\)+1;end) s
( ) s
( for ifreq=ilow:ihigh) s
( ) s
( f = 2.*pi*\(ifreq-1\)/nt/dt;) s
( L = exp\(i*f*\(h.^N\)'*q\); ) s
( y = D\(ifreq,:\)';) s
( x = L'*y;) s
( ) s
( MATRIX = L'*L;) s
(  tr=real\(trace\(MATRIX\)\);) s
(  Q =mu*tr*eye\(nq\);) s
(  x = inv\(MATRIX+Q\) *L'* y; ) s
( M\(ifreq,:\) = x';) s
( M\(nt+2-ifreq,:\) = conj\(x\)';) s
( end) s
( ) s
( M\(nt/2+1,:\) = zeros\(1,nq\);) s
( m = real\(ifft\(M,[],1\)\);) s
( return) s
( ) s
/sd 1 def
/sd 0 def
1 sn
fnf
lg lgp
pagesave restore
showpage

%%Trailer
%%Pages: 1
docsave restore end