Home > atmlab > demos > asg2y_demo2.m

asg2y_demo2

PURPOSE ^

ASG2Y_DEMO2 A simple demonstration how to use

SYNOPSIS ^

function [Y,P]=asg2y_demo2(sensor,G,FS)

DESCRIPTION ^

 ASG2Y_DEMO2   A simple demonstration how to use
    asg functions to simulate radiative transfer
    for generated atmospheric states (by asg_demo)

    Gives an example on basic usage of ASG2Y

    The example simulate radiative transfer 
    for a cloudy state generated by asg_demo 

    A lot of default settings are made in this
    demo, which may be re-considered depending
    on what data should be used for

 FORMAT [Y,P]=asg2y_demo2(sensor,G,FS)

 OUT    
         Y      brightness temperature (planck Tb)
         P      antenna weighted atmospheric profiles       

 IN      sensor structure with fields
                name: i.e amsu or mhs
                channels: [1-5]
                los: center line of sights
                pos: sensor position 
         G      structure with atmospheric variables on gformat data
                as produced by asg_demo2
         FS     FS   structure specifying atmospheric settings      
                as produced by asg_demo2

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

asg2y_demo2.m

SOURCE CODE ^

0001 % ASG2Y_DEMO2   A simple demonstration how to use
0002 %    asg functions to simulate radiative transfer
0003 %    for generated atmospheric states (by asg_demo)
0004 %
0005 %    Gives an example on basic usage of ASG2Y
0006 %
0007 %    The example simulate radiative transfer
0008 %    for a cloudy state generated by asg_demo
0009 %
0010 %    A lot of default settings are made in this
0011 %    demo, which may be re-considered depending
0012 %    on what data should be used for
0013 %
0014 % FORMAT [Y,P]=asg2y_demo2(sensor,G,FS)
0015 %
0016 % OUT
0017 %         Y      brightness temperature (planck Tb)
0018 %         P      antenna weighted atmospheric profiles
0019 %
0020 % IN      sensor structure with fields
0021 %                name: i.e amsu or mhs
0022 %                channels: [1-5]
0023 %                los: center line of sights
0024 %                pos: sensor position
0025 %         G      structure with atmospheric variables on gformat data
0026 %                as produced by asg_demo2
0027 %         FS     FS   structure specifying atmospheric settings
0028 %                as produced by asg_demo2
0029 
0030 % Created by Bengt Rydberg
0031 
0032 
0033 function [Y,P]=asg2y_demo2(sensor,G,FS)
0034 
0035   error( 'Can not be used presently. Not yet updated.' );
0036 
0037 atmlab('FMODEL_VERBOSITY',0)
0038 
0039 %-----------------------------------
0040 %SET Q-STUCTURE (for ARTS RT-simulations)
0041 %in sub-function asg_q a lot of default settings are made
0042 [Q,sensor,FS]=asg_q(sensor,FS);
0043 
0044 %------------------------------------
0045 %ADD SCATTERING PROPERTIES OF ICE PARTICLES FOR THE RELEVANT
0046 %INSTRUMENT FREQUENCIES
0047 G=asg_ssp(Q.F_GRID,G);
0048       
0049 %---------------------------------------
0050 %RUN RTE CALCULATIONS
0051 workfolder =create_tmpfolder;
0052 [Y]=asg_gfs2i(G,Q,FS,sensor,workfolder);
0053 delete_tmpfolder(workfolder)
0054 
0055 %--------------------------------------
0056 %extract atmospheric parameters in line of sights
0057 P=asg_atm_pb(G,Q,sensor)
0058 
0059 
0060 %-sub-functions
0061 
0062 
0063 function [Q,instrument,FS]=asg_q(instrument,FS)
0064 
0065  Q = qarts;
0066  Q.OUTPUT_FILE_FORMAT='binary';
0067 
0068  %- Basic atmospheric variables
0069  %
0070  Q.ATMOSPHERE_DIM = 3;
0071  Q.P_GRID         = 100/16e3; % Unit here is pressure decade
0072  Q.LAT_GRID=[];
0073  Q.LON_GRID=[];
0074  Q.R_GEOID=6366700;
0075 
0076 
0077 
0078  %- Spectroscopic stuff
0079  %
0080  Q.STOKES_DIM       = 1;
0081  Q.ABS_LINES_FORMAT = 'Arts';
0082  Q.ABSORPTION='OnTheFly';
0083  Q.ABS_LINESHAPE='Voigt_Kuntz6';
0084  Q.ABS_LINESHAPE_CUTOFF=-1;
0085  Q.ABS_LINESHAPE_FACTOR='quadratic';
0086  
0087  %- RTE variables
0088  %
0089  Q.PPATH_STEP_AGENDA = { 'ppath_stepGeometric( ppath_step, atmosphere_dim, p_grid, lat_grid, lon_grid, z_field, r_geoid, z_surface, ppath_lmax)' };
0090  Q.PPATH_LMAX=1e3;
0091  Q.RTE_AGENDA= {'RteStd(iy,diy_dvmr,diy_dt,ppath,ppath_array,ppath_array_index,f_grid,stokes_dim,emission_agenda,abs_scalar_gas_agenda,rte_do_vmr_jacs,rte_do_t_jacs )'};
0092  Q.EMISSION_AGENDA={'emissionPlanck'};
0093  Q.IY_SPACE_AGENDA={'Ignore(rte_los)','Ignore(rte_pos)',...
0094               'MatrixCBR(iy,stokes_dim,f_grid)' };
0095  Q.J_DO=0;
0096 
0097  %- SENSOR variabels
0098  Q.ANTENNA_DIM=1;
0099  Q.IY_UNIT            = 'RJBT';
0100  Q.SENSOR_DO=0;
0101 
0102  %Sensor stuff
0103 
0104   
0105   sensor_los=instrument.los;
0106   sensor_pos=instrument.pos;
0107   m_za_grid=[-1:0.25:1];
0108   m_az_grid=[-1:0.25:1];   
0109   for i=1:length(m_az_grid)
0110        [za_c(:,i),aa_c(:,i)] = ...
0111        arts_map_daa( sensor_los(1)+m_za_grid(i), sensor_los(2), m_az_grid );
0112   end
0113       
0114   Q.SENSOR_LOS=zeros(length(za_c(:)),2);
0115   Q.SENSOR_LOS(:,1)=za_c(:);
0116   Q.SENSOR_LOS(:,2)=aa_c(:);
0117   Q.SENSOR_POS=zeros(length(za_c(:)),3);
0118   Q.SENSOR_POS(:,1)=sensor_pos(1);
0119   Q.SENSOR_POS(:,2)=sensor_pos(2);
0120   Q.SENSOR_POS(:,3)=sensor_pos(3);  
0121    
0122  if strcmp('amsu',instrument.name)
0123 
0124   f_grid_file=fullfile( atmlab('ARTS_INCLUDES'),'amsu/amsub.f_grid_fast.xml');
0125   sensor_response_file= ...
0126   fullfile( atmlab('ARTS_INCLUDES'),'amsu/amsub.sensor_response_fast.xml');
0127  
0128  end
0129  if strcmp('mhs',instrument.name)
0130 
0131   f_grid_file=fullfile( atmlab('ARTS_INCLUDES'),'mhs/mhs.f_grid_fast.xml');
0132   sensor_response_file= ...
0133   fullfile( atmlab('ARTS_INCLUDES'),'mhs/mhs.sensor_response_fast.xml');
0134     
0135  end
0136 
0137   f_grid=xmlLoad(f_grid_file);
0138 
0139   s_r=full(xmlLoad(sensor_response_file));
0140  
0141   H = amsu_antenna(m_za_grid,m_az_grid);
0142   R=full(H);
0143   ind=[];
0144   k=1;
0145   for i=instrument.channels
0146       ind=[ind,find(s_r(i,:))];
0147   end
0148   ind1=sort(ind);
0149   instrument.ch_resp=s_r(instrument.channels,ind1);
0150 
0151   ind2=[];
0152   for j=1:length(ind)
0153       ind2=[ind2,[ind1(j):length(f_grid):length(R)]];
0154   end
0155   ind2=sort(ind2);
0156   instrument.ant_resp=R(ind1,ind2);
0157 
0158   Q.F_GRID           = f_grid(ind1);
0159   if strcmp('amsu',instrument.name) | strcmp('amsu',instrument.name)
0160      Q.ABS_LINES=...
0161      fullfile(atmlab('ARTS_INCLUDES'),'amsu/amsub.hitran04_lines.xml');
0162   end
0163   if strcmp('mhs',instrument.name)
0164      Q.ABS_LINES=...
0165      fullfile(atmlab('ARTS_INCLUDES'),'mhs/mhs.hitran04_lines.xml');
0166   end
0167  
0168 
0169   Q.ABS_MODELS={fullfile(atmlab('ARTS_INCLUDES'),'continua.arts')};
0170  
0171 
0172  %- Surface properties
0173  %
0174  if strcmp('amsu',instrument.name) | strcmp('mhs',instrument.name)
0175 
0176    %randomize surface emissivity
0177    emis=min(0.95+0.05*randn(length(Q.F_GRID),1),1);
0178    emis=max(0.8,emis);
0179    FS.SURF_EMIS=emis;
0180    %SURFACE EMISSIVITIES
0181    Se=[];
0182    for j=1:length(FS.SURF_EMIS)
0183     se=num2str(FS.SURF_EMIS(j));
0184     Se=[Se,',',se];
0185    end  
0186    Se=Se(2:end);     
0187    Q.SURFACE_PROP_AGENDA= ...
0188     {['InterpAtmFieldToRteGps(surface_skin_t,atmosphere_dim,'...
0189       'p_grid, lat_grid, lon_grid,rte_gp_p, rte_gp_lat,'...
0190       'rte_gp_lon,t_field)'],...
0191       'Ignore(rte_pos)',...
0192       ['VectorCreate(surface_emissivity)'],... 
0193      ['VectorSet(surface_emissivity,[',Se,'])'],...
0194       ['surfaceFlatVaryingEmissivity(surface_los, surface_rmatrix,'...
0195       'surface_emission, f_grid, stokes_dim,atmosphere_dim, rte_los,'...
0196       'surface_skin_t,surface_emissivity)']}; 
0197       
0198   end
0199 
0200  %- MC variables
0201  %
0202  %antenna response for one channel, will be used to setup
0203  %optimized monte carlo noise
0204  R=R(ind1(1),ind1(1):length(s_r):end); 
0205  std1=0.5;
0206  std_max=6;
0207  [mc_std]=mc_error(R,std1,std_max);
0208 
0209  Q.CLOUDBOX.METHOD                      = 'MC';
0210  Q.CLOUDBOX.METHOD_PRMTRS.STD_ERR       = mc_std;      
0211  Q.CLOUDBOX.METHOD_PRMTRS.MAX_TIME      = -1;
0212  Q.CLOUDBOX.METHOD_PRMTRS.MAX_ITER      = -1;
0213  Q.CLOUDBOX.METHOD_PRMTRS.Z_FIELD_IS_1D = 0;
0214  Q.CLOUDBOX.LIMITS                      = [];
0215  Q.CLOUDBOX.OPT_PROP_GAS_AGENDA={'ext_matInit','abs_vecInit',...
0216                                 'ext_matAddGas','abs_vecAddGas'};
0217  Q.CLOUDBOX.OPT_PROP_PART_AGENDA={'ext_matAddPart','abs_vecAddPart'};
0218  Q.CLOUDBOX.SPT_CALC_AGENDA={};
0219  Q.CLOUDBOX_DO=0;
0220 
0221 return
0222 
0223 %-------------------------------------------------------------------------
0224 
0225 function G=asg_ssp(f_grid,G)
0226  %properties for adding pnd fields
0227  alpha=0;Ni=10;xnorm=50e-6;
0228  [x_i,w_i]=gauss_laguerre(alpha,Ni,xnorm);
0229  Nj=3;xnormw=10e-6;
0230  [x_k,w_k]=gauss_laguerre(alpha,Nj,xnormw);
0231  %x_i will be the diameter of the particles considered
0232  %calculate ssp
0233  for i=1:2
0234 
0235   if i==1 
0236    F=create_ssp('sphere',f_grid,x_i/2);
0237    P(i).x=x_i;
0238    P(i).w=w_i;
0239    P(i).x_norm=xnorm;
0240    P(i).PSD='MH97';
0241   else
0242    T_grid=[273 283 293 303 313];
0243    for i1=1:length(f_grid)
0244      for j1=1:length(T_grid)
0245        rfr_index(i1,j1) = sqrt(eps_water_liebe93( f_grid(i1), T_grid(j1) ));
0246      end
0247    end
0248    F=create_ssp('sphere',f_grid,x_k,T_grid,rfr_index);
0249    P(i).x=x_k;
0250    P(i).w=w_k;
0251    P(i).x_norm=xnormw;
0252    P(i).PSD='Water';
0253   end
0254   P(i).SSP=F;
0255   P(i).method='gauss-laguerre';
0256   P(i).shape='sphere';  
0257  end
0258 
0259  ind=find(strncmp({G.NAME},'Particles',9));
0260  for j=ind
0261      G(j).PROPS=P(1).SSP(j-ind(1)+1);
0262  end
0263 
0264 
0265 return
0266 
0267 
0268 %-------------------------------------------------------
0269 function [Y]=asg_gfs2i(G,Q,FS,sensor,workfolder)
0270  
0271  Q.P_GRID=G(1).GRID1;
0272  Q.LAT_GRID=G(1).GRID2;
0273  Q.LON_GRID=G(1).GRID3;
0274 
0275 
0276  %- Determine necessary size of cloud box and adjust it
0277  %
0278  [Q,G]=asg_set_cb(Q,G);
0279 
0280 
0281  %check if there is any liquid cloud
0282  %if not we don't need to include it in the simulations
0283  ind=find(strcmp({G.NAME},'LWC field'));
0284  if ~isempty(ind)
0285    ind=find(strcmp({G.NAME},'LWC field'));
0286    G(ind).PROPS=['liquidcloud-MPM93'];
0287    G(ind).DATA=G(ind).DATA*1e-3;
0288    G(ind).DATA_NAME='Volume mixing ratio';
0289    G(ind)= asg_regrid( asgD, G(ind), Q );
0290    G(ind).DATA(find(G(ind).DATA(:)>10e-3))=10e-3;
0291  else
0292    ind=setdiff(1:length(G),ind);
0293    G=G(ind);
0294  end
0295 
0296  %- Run ARTS
0297  %
0298  if strcmp(sensor.name,'amsu') | strcmp(sensor.name,'mhs')
0299   
0300    %create an include file
0301    fid=fopen(fullfile(workfolder,'include.arts'));
0302    if fid==-1
0303       fid=fopen(fullfile(workfolder,'include.arts'),'a');
0304       fprintf(fid,'Arts2 {\n'); 
0305       fprintf(fid,'#VectorSetExplicitly(abs_n2,[0.7808])\n');
0306       fprintf(fid,'VectorSet( abs_n2, [ 0.7808 ] )\n');
0307       fprintf(fid,'IndexSet( abs_p_interp_order,   5 )\n');
0308       fprintf(fid,'IndexSet( abs_t_interp_order,   7 )\n');
0309       fprintf(fid,'IndexSet( abs_nls_interp_order, 5 )}\n');
0310       fclose(fid);
0311    end
0312    Q.INCLUDES= {fullfile(workfolder,'include.arts')};
0313    ind=find(Q.SENSOR_LOS(:,2)>180);
0314    if length(ind)>0
0315       Q.SENSOR_LOS(ind,2)=Q.SENSOR_LOS(ind,2)-360;
0316    end 
0317      
0318    if Q.CLOUDBOX_DO==0 %use clearsky mode
0319 
0320       Q = asg2q( asgD, G, Q, workfolder );  
0321       Q.ABS_NLS=[];
0322       Q                  = qarts_abstable( Q,9,100,[],[] );
0323       [y1,ydata,dy] = arts_y( Q, workfolder );
0324 
0325    else %do MC simulations
0326      
0327      y1=zeros(length(Q.F_GRID)*length(Q.SENSOR_LOS),1);
0328      std_err=unique(Q.CLOUDBOX.METHOD_PRMTRS.STD_ERR);
0329    
0330      %loop over different part of the "antenna"
0331      Q1=Q; 
0332      for i=1:length(std_err)
0333       i
0334       ind=find(Q.CLOUDBOX.METHOD_PRMTRS.STD_ERR==std_err(i));
0335       Q1.CLOUDBOX.METHOD_PRMTRS.STD_ERR=std_err(i);
0336       if i==1
0337        Q1 = asg2q( asgD, G, Q1, workfolder );  
0338        Q1.ABS_NLS=[];
0339        Q1 = qarts_abstable( Q1,9,100,[],[] );
0340       end
0341 
0342       Q1.SENSOR_POS=Q.SENSOR_POS(ind,:);
0343       Q1.SENSOR_LOS=Q.SENSOR_LOS(ind,:);
0344     
0345       [y,ydata,dy] = arts_y( Q1, workfolder );
0346     
0347       for j=1:length(ind)
0348         y1((ind(j)-1)*length(Q.F_GRID)+1:ind(j)*length(Q.F_GRID))=...
0349         y((j-1)*length(Q.F_GRID)+1:j*length(Q.F_GRID));
0350       end 
0351    
0352      end
0353  
0354    end 
0355  end
0356    
0357  %convert to intensity
0358  boltzmann = constants('BOLTZMANN_CONST');
0359  speed_light = constants('SPEED_OF_LIGHT');
0360  f_grid=zeros(size(y1));
0361  for j=1:length(Q.F_GRID)
0362     f_grid(j:length(Q.F_GRID):end)=Q.F_GRID(j);
0363  end
0364  %convert to Planck
0365    
0366  I=2*f_grid.^2*boltzmann/(speed_light^2).*y1;
0367  Y= i2planckTb(I,f_grid);
0368 
0369  %weight pencil beam simulations with antenna and sideband response
0370  Y=sensor.ch_resp*sensor.ant_resp*Y;
0371 
0372  
0373 return
0374 
0375 %-----------------------------------------------------------------
0376 
0377 function [Q,G]=asg_set_cb(Q,G)
0378 
0379 iwc_ind = find( strcmp( lower({G.NAME}), 'iwc field' ) );
0380 [C1,C2,C3,C4,C5,C6]=deal([]);
0381 
0382 for i=1:size(G(iwc_ind).DATA,3)
0383   [c1,c2]=find(G(iwc_ind).DATA(:,:,i));
0384   if ~isempty(c1)
0385      C1=min([C1,vec2col(c1)']);
0386      C2=max([C2,vec2col(c1)']);
0387   end
0388   if ~isempty(c2)
0389      C3=min([C3,vec2col(c2)']);
0390      C4=max([C4,vec2col(c2)']);
0391   end
0392 end
0393 
0394 
0395 for i=1:size(G(iwc_ind).DATA,1)
0396   [c1,c2]=find(squeeze(G(iwc_ind).DATA(i,:,:)));
0397     if ~isempty(c2)
0398        C5=min([C5,vec2col(c2)']);
0399        C6=max([C6,vec2col(c2)']);
0400    end
0401 end
0402 %
0403 
0404 if ~isempty(C1) %if C1 is empty we don't have any cloud particles
0405     p1   = G(iwc_ind).GRID1(C1);
0406     p2   = G(iwc_ind).GRID1(C2);
0407     lat1 = G(iwc_ind).GRID2(C3);
0408     lat2 = G(iwc_ind).GRID2(C4);
0409     z    = p2z_cira86([p1 p2]',mean([lat1 lat2]),160);
0410     lon1 = G(iwc_ind).GRID3(C5);
0411     lon2 = G(iwc_ind).GRID3(C6);
0412     %
0413     Q.CLOUDBOX.LIMITS = [z(1)-2e3 z(2)+2e3 lat1 lat2 lon1 lon2];
0414     Q.CLOUDBOX_DO=1;
0415 else
0416     %remove all particles in G
0417     ind  =find(strncmp({G.NAME},'Particles',9)==0);
0418     G    = G(ind);
0419 
0420 end
0421 
0422 return
0423 
0424 %----------------------------------------
0425 function P=asg_atm_pb(G,Q,sensor)
0426  
0427  addpath /home/bengt/ARTS/atmlab/physics/thermodynamic
0428 
0429  p = asg_pathiwc( asgD, G, Q, 21e3 );
0430  zvec=[0:250:18e3];
0431  warning off
0432  for i=1:length(p)
0433     iwc(:,i)=interp1(p(i).z,p(i).iwc,zvec);
0434     lwc(:,i)=interp1(p(i).z,p(i).lwc,zvec);
0435     iwp(i)=p(i).iwp;
0436     rhi(:,i)=interp1(p(i).z,p(i).rhi,zvec);
0437     t(:,i)=interp1(p(i).z,p(i).t,zvec);
0438  end
0439  warning on
0440  %weight the individual profiles
0441  w=sensor.ant_resp(1,1:length(Q.F_GRID):end);
0442  P.z=zvec;
0443  P.iwc=iwc*w';
0444  P.lwc=lwc*w';
0445  P.iwp=iwp*w';
0446  P.rhi=rhi*w';
0447  P.t=t*w';
0448 
0449 return
0450 
0451 function H = amsu_antenna(za_grid,aa_grid)
0452 
0453  Q = qarts;
0454 
0455 
0456  Q.SENSOR_DO      = 1;
0457  Q.STOKES_DIM     = 1;
0458  Q.ANTENNA_DIM    = 2;
0459  Q.ATMOSPHERE_DIM = 3; 
0460 
0461  Q.OUTPUT_FILE_FORMAT = 'ascii';
0462 
0463  Q.F_GRID = ...
0464     fullfile( atmlab('ARTS_INCLUDES'), 'amsu', 'amsub.f_grid_fast.xml' );
0465 
0466  Q.MBLOCK_ZA_GRID = za_grid;
0467  Q.MBLOCK_AA_GRID = aa_grid;
0468 
0469  Q.SENSOR_RESPONSE = qartsSensor;
0470 
0471  Q.SENSOR_RESPONSE.SENSOR_NORM      = 1;
0472  Q.SENSOR_RESPONSE.ANTENNA_LOS      = [0 0];
0473  Q.SENSOR_RESPONSE.ANTENNA_RESPONSE = fullfile( atmlab('ARTS_INCLUDES'), ...
0474                                                'amsu', 'amsub.antenna.xml' );
0475  Q = arts_sensor( Q );
0476 
0477  H = Q.SENSOR_RESPONSE;
0478 
0479 return
0480 
0481 %-------
0482 
0483 function [mc_std]=mc_error(w,std,std_max)
0484 
0485 
0486  mc_std=sqrt(std^2./w);
0487 
0488  mc_std(find(mc_std>std_max))=std_max;
0489 
0490  ind1=find(mc_std(:)==std_max); 
0491  ind1b=ind1;
0492  ind2=find(mc_std(:)<std_max); 
0493 
0494  while length(ind1)>0
0495   varp=std^2-sum(w(ind1b).^2.*mc_std(ind1b).^2);
0496   mc_std(ind2)=min(sqrt(varp/sum(w(ind2))./w(ind2)),std_max);
0497   ind1=find(mc_std(ind2)>std_max);
0498   if length(ind1)>0
0499      ind1b=sort([ind1b,ind2(ind1)])
0500   end
0501   ind2=setdiff(1:length(mc_std(:)),ind1b);   
0502  end
0503 
0504 return

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