Home > atmlab > forwardmodel > voigt_humlik.m

voigt_humlik

PURPOSE ^

VOIGT_HUMLIK Voigt line shape

SYNOPSIS ^

function bredd=voigt_humlik(f,f0,dfd,dfp)

DESCRIPTION ^

 VOIGT_HUMLIK   Voigt line shape 

    A simple Matlab implementation of the Humlicek algorithm for
    calculation the Voigt line shape (JQSRT, 21, 309-313, 1978).

    Please note that width for both Doppler and pressure broadening is given
    as FWHM/2. This corresponds to the standard width for pressure
    broadening, but could differ for the Doppler part. See *doppler_width*
    for some details.

    The function handles only single altitudes (but multiple frequencies).

 FORMAT   bredd = voigt_humlik(f,f0,dfd,dfp)
        
 OUT   bredd   Line shape.
 IN    f       Frequency vector.
       f0      Centre frequency
       dfd     Width of Doppler broadening.
       dfp     Width of pressure broadening.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

DOWNLOAD ^

voigt_humlik.m

SOURCE CODE ^

0001 % VOIGT_HUMLIK   Voigt line shape
0002 %
0003 %    A simple Matlab implementation of the Humlicek algorithm for
0004 %    calculation the Voigt line shape (JQSRT, 21, 309-313, 1978).
0005 %
0006 %    Please note that width for both Doppler and pressure broadening is given
0007 %    as FWHM/2. This corresponds to the standard width for pressure
0008 %    broadening, but could differ for the Doppler part. See *doppler_width*
0009 %    for some details.
0010 %
0011 %    The function handles only single altitudes (but multiple frequencies).
0012 %
0013 % FORMAT   bredd = voigt_humlik(f,f0,dfd,dfp)
0014 %
0015 % OUT   bredd   Line shape.
0016 % IN    f       Frequency vector.
0017 %       f0      Centre frequency
0018 %       dfd     Width of Doppler broadening.
0019 %       dfp     Width of pressure broadening.
0020 
0021 % 1992         Mats Pettersson
0022 % 1993         Modified by Magnus Gustafsson
0023 % 2006-11-26   Header written by Patrick Eriksson.
0024 
0025 
0026 function bredd=voigt_humlik(f,f0,dfd,dfp)
0027 
0028 bl=dfp;        %% MG
0029 bd=dfd;        %% MG
0030 x=f-f0;        %% MG
0031 
0032 
0033 
0034 T=[.314240376 .947788391 1.59768264 2.27950708 3.02063703 3.8897249];
0035 C=[1.01172805 -.75197147 1.2557727e-2 1.00220082e-2 -2.42068135e-4 5.00848061e-7];
0036 S=[1.393237 .231152406 -.155351466 6.21836624e-3 9.19082986e-5 -6.27525958e-7];
0037 
0038 
0039 p=0;
0040 WR=0;WR2=0;WR1=0;
0041 WI=0;
0042 X=x/bd*(0.83255461115770);
0043 Y=bl/bd*(0.83255461115770);
0044 
0045 Y1=Y+1.5;
0046 Y2=Y1*Y1;
0047 X1=X;
0048 X2=X;
0049 n=0;
0050 m=0;
0051 k=0;
0052 l=0;
0053 if Y>0.85
0054 
0055  for i=1:6
0056    R=X-T(i);
0057    D=(R.^2+Y2).^(-1);
0058    D1=Y1*D;
0059    D2=R.*D;
0060    R=X+T(i);
0061    D=(R.^2+Y2).^(-1);
0062    D3=Y1*D;
0063    D4=R.*D;
0064    WR=WR+C(i)*(D1+D3)-S(i)*(D2-D4);
0065    
0066  end
0067 
0068 else
0069 
0070  vektorlangd=length(X);
0071  test=sign(abs(X)-(18.1*Y+1.65));
0072  
0073  for p=1:vektorlangd
0074    if test(p)<0
0075      k=k+1;
0076      X1(k)=X(p);
0077    else
0078      n=n+1;
0079      X2(n)=X(p);
0080    end
0081  end
0082  X1=X1(1:k);
0083  X2=X2(1:n);
0084 
0085  if n>0
0086    WR2=exp(-X2.*X2);
0087    Y3=Y+3;
0088    for i=1:6
0089      R=X2-T(i);
0090      D=(R.^2+Y2).^(-1);
0091      D1=Y1*D;
0092      D2=R.*D;
0093      WR2=WR2+Y*(C(i)*(R.*D2-1.5*D1)+S(i)*Y3*D2).*((R.^2+2.25).^(-1));
0094      R=X2+T(i);
0095      D=(R.^2+Y2).^(-1);
0096      D3=Y1*D;
0097      D4=R.*D;
0098      WR2=WR2+Y*(C(i)*(R.*D4-1.5*D3)-S(i)*Y3*D4).*((R.^2+2.25).^(-1));  
0099     
0100    end
0101  
0102  end
0103 
0104  if k>0
0105    for i=1:6
0106     R=X1-T(i);
0107     D=(R.^2+Y2).^(-1);
0108     D1=Y1*D;
0109     D2=R.*D;
0110     R=X1+T(i);
0111     D=(R.^2+Y2).^(-1);
0112     D3=Y1*D;
0113     D4=R.*D;
0114     WR1=WR1+C(i)*(D1+D3)-S(i)*(D2-D4);
0115     
0116    end
0117  end  
0118  
0119  for p=1:vektorlangd
0120    if test(p)<0
0121      l=l+1;
0122      WR(p)=WR1(l);
0123    else
0124      m=m+1;
0125      WR(p)=WR2(m);
0126    end
0127  end
0128 
0129 end  
0130     
0131 
0132 %----Modification !-----
0133 s=size(WR);
0134 if s(1)<s(2)
0135     WR=WR';
0136 end
0137 %-----------------------
0138  bredd=1/bd*(0.46971863934983)*WR;
0139 
0140 
0141 end

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005