% 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 % Created by Bengt Rydberg function [Y,P]=asg2y_demo2(sensor,G,FS) error( 'Can not be used presently. Not yet updated.' ); atmlab('FMODEL_VERBOSITY',0) %----------------------------------- %SET Q-STUCTURE (for ARTS RT-simulations) %in sub-function asg_q a lot of default settings are made [Q,sensor,FS]=asg_q(sensor,FS); %------------------------------------ %ADD SCATTERING PROPERTIES OF ICE PARTICLES FOR THE RELEVANT %INSTRUMENT FREQUENCIES G=asg_ssp(Q.F_GRID,G); %--------------------------------------- %RUN RTE CALCULATIONS workfolder =create_tmpfolder; [Y]=asg_gfs2i(G,Q,FS,sensor,workfolder); delete_tmpfolder(workfolder) %-------------------------------------- %extract atmospheric parameters in line of sights P=asg_atm_pb(G,Q,sensor) %-sub-functions function [Q,instrument,FS]=asg_q(instrument,FS) Q = qarts; Q.OUTPUT_FILE_FORMAT='binary'; %- Basic atmospheric variables % Q.ATMOSPHERE_DIM = 3; Q.P_GRID = 100/16e3; % Unit here is pressure decade Q.LAT_GRID=[]; Q.LON_GRID=[]; Q.R_GEOID=6366700; %- Spectroscopic stuff % Q.STOKES_DIM = 1; Q.ABS_LINES_FORMAT = 'Arts'; Q.ABSORPTION='OnTheFly'; Q.ABS_LINESHAPE='Voigt_Kuntz6'; Q.ABS_LINESHAPE_CUTOFF=-1; Q.ABS_LINESHAPE_FACTOR='quadratic'; %- RTE variables % Q.PPATH_STEP_AGENDA = { 'ppath_stepGeometric( ppath_step, atmosphere_dim, p_grid, lat_grid, lon_grid, z_field, r_geoid, z_surface, ppath_lmax)' }; Q.PPATH_LMAX=1e3; 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 )'}; Q.EMISSION_AGENDA={'emissionPlanck'}; Q.IY_SPACE_AGENDA={'Ignore(rte_los)','Ignore(rte_pos)',... 'MatrixCBR(iy,stokes_dim,f_grid)' }; Q.J_DO=0; %- SENSOR variabels Q.ANTENNA_DIM=1; Q.Y_UNIT = 'RJBT'; Q.SENSOR_DO=0; %Sensor stuff sensor_los=instrument.los; sensor_pos=instrument.pos; m_za_grid=[-1:0.25:1]; m_az_grid=[-1:0.25:1]; for i=1:length(m_az_grid) [za_c(:,i),aa_c(:,i)] = ... arts_map_daa( sensor_los(1)+m_za_grid(i), sensor_los(2), m_az_grid ); end Q.SENSOR_LOS=zeros(length(za_c(:)),2); Q.SENSOR_LOS(:,1)=za_c(:); Q.SENSOR_LOS(:,2)=aa_c(:); Q.SENSOR_POS=zeros(length(za_c(:)),3); Q.SENSOR_POS(:,1)=sensor_pos(1); Q.SENSOR_POS(:,2)=sensor_pos(2); Q.SENSOR_POS(:,3)=sensor_pos(3); if strcmp('amsu',instrument.name) f_grid_file=fullfile( atmlab('ARTS_INCLUDES'),'amsu/amsub.f_grid_fast.xml'); sensor_response_file= ... fullfile( atmlab('ARTS_INCLUDES'),'amsu/amsub.sensor_response_fast.xml'); end if strcmp('mhs',instrument.name) f_grid_file=fullfile( atmlab('ARTS_INCLUDES'),'mhs/mhs.f_grid_fast.xml'); sensor_response_file= ... fullfile( atmlab('ARTS_INCLUDES'),'mhs/mhs.sensor_response_fast.xml'); end f_grid=xmlLoad(f_grid_file); s_r=full(xmlLoad(sensor_response_file)); H = amsu_antenna(m_za_grid,m_az_grid); R=full(H); ind=[]; k=1; for i=instrument.channels ind=[ind,find(s_r(i,:))]; end ind1=sort(ind); instrument.ch_resp=s_r(instrument.channels,ind1); ind2=[]; for j=1:length(ind) ind2=[ind2,[ind1(j):length(f_grid):length(R)]]; end ind2=sort(ind2); instrument.ant_resp=R(ind1,ind2); Q.F_GRID = f_grid(ind1); if strcmp('amsu',instrument.name) | strcmp('amsu',instrument.name) Q.ABS_LINES=... fullfile(atmlab('ARTS_INCLUDES'),'amsu/amsub.hitran04_lines.xml'); end if strcmp('mhs',instrument.name) Q.ABS_LINES=... fullfile(atmlab('ARTS_INCLUDES'),'mhs/mhs.hitran04_lines.xml'); end Q.ABS_MODELS={fullfile(atmlab('ARTS_INCLUDES'),'continua.arts')}; %- Surface properties % if strcmp('amsu',instrument.name) | strcmp('mhs',instrument.name) %randomize surface emissivity emis=min(0.95+0.05*randn(length(Q.F_GRID),1),1); emis=max(0.8,emis); FS.SURF_EMIS=emis; %SURFACE EMISSIVITIES Se=[]; for j=1:length(FS.SURF_EMIS) se=num2str(FS.SURF_EMIS(j)); Se=[Se,',',se]; end Se=Se(2:end); Q.SURFACE_PROP_AGENDA= ... {['InterpAtmFieldToRteGps(surface_skin_t,atmosphere_dim,'... 'p_grid, lat_grid, lon_grid,rte_gp_p, rte_gp_lat,'... 'rte_gp_lon,t_field)'],... 'Ignore(rte_pos)',... ['VectorCreate(surface_emissivity)'],... ['VectorSet(surface_emissivity,[',Se,'])'],... ['surfaceFlatVaryingEmissivity(surface_los, surface_rmatrix,'... 'surface_emission, f_grid, stokes_dim,atmosphere_dim, rte_los,'... 'surface_skin_t,surface_emissivity)']}; end %- MC variables % %antenna response for one channel, will be used to setup %optimized monte carlo noise R=R(ind1(1),ind1(1):length(s_r):end); std1=0.5; std_max=6; [mc_std]=mc_error(R,std1,std_max); Q.CLOUDBOX.METHOD = 'MC'; Q.CLOUDBOX.METHOD_PRMTRS.STD_ERR = mc_std; Q.CLOUDBOX.METHOD_PRMTRS.MAX_TIME = -1; Q.CLOUDBOX.METHOD_PRMTRS.MAX_ITER = -1; Q.CLOUDBOX.METHOD_PRMTRS.Z_FIELD_IS_1D = 0; Q.CLOUDBOX.LIMITS = []; Q.CLOUDBOX.OPT_PROP_GAS_AGENDA={'ext_matInit','abs_vecInit',... 'ext_matAddGas','abs_vecAddGas'}; Q.CLOUDBOX.OPT_PROP_PART_AGENDA={'ext_matAddPart','abs_vecAddPart'}; Q.CLOUDBOX.SPT_CALC_AGENDA={}; Q.CLOUDBOX_DO=0; return %------------------------------------------------------------------------- function G=asg_ssp(f_grid,G) %properties for adding pnd fields alpha=0;Ni=10;xnorm=50e-6; [x_i,w_i]=gauss_laguerre(alpha,Ni,xnorm); Nj=3;xnormw=10e-6; [x_k,w_k]=gauss_laguerre(alpha,Nj,xnormw); %x_i will be the diameter of the particles considered %calculate ssp for i=1:2 if i==1 F=create_ssp('sphere',f_grid,x_i/2); P(i).x=x_i; P(i).w=w_i; P(i).x_norm=xnorm; P(i).PSD='MH97'; else T_grid=[273 283 293 303 313]; for i1=1:length(f_grid) for j1=1:length(T_grid) rfr_index(i1,j1) = sqrt(eps_water_liebe93( f_grid(i1), T_grid(j1) )); end end F=create_ssp('sphere',f_grid,x_k,T_grid,rfr_index); P(i).x=x_k; P(i).w=w_k; P(i).x_norm=xnormw; P(i).PSD='Water'; end P(i).SSP=F; P(i).method='gauss-laguerre'; P(i).shape='sphere'; end ind=find(strncmp({G.NAME},'Particles',9)); for j=ind G(j).PROPS=P(1).SSP(j-ind(1)+1); end return %------------------------------------------------------- function [Y]=asg_gfs2i(G,Q,FS,sensor,workfolder) Q.P_GRID=G(1).GRID1; Q.LAT_GRID=G(1).GRID2; Q.LON_GRID=G(1).GRID3; %- Determine necessary size of cloud box and adjust it % [Q,G]=asg_set_cb(Q,G); %check if there is any liquid cloud %if not we don't need to include it in the simulations ind=find(strcmp({G.NAME},'LWC field')); if ~isempty(ind) ind=find(strcmp({G.NAME},'LWC field')); G(ind).PROPS=['liquidcloud-MPM93']; G(ind).DATA=G(ind).DATA*1e-3; G(ind).DATA_NAME='Volume mixing ratio'; G(ind)= asg_regrid( asgD, G(ind), Q ); G(ind).DATA(find(G(ind).DATA(:)>10e-3))=10e-3; else ind=setdiff(1:length(G),ind); G=G(ind); end %- Run ARTS % if strcmp(sensor.name,'amsu') | strcmp(sensor.name,'mhs') %create an include file fid=fopen(fullfile(workfolder,'include.arts')); if fid==-1 fid=fopen(fullfile(workfolder,'include.arts'),'a'); fprintf(fid,'Arts2 {\n'); fprintf(fid,'#VectorSetExplicitly(abs_n2,[0.7808])\n'); fprintf(fid,'VectorSet( abs_n2, [ 0.7808 ] )\n'); fprintf(fid,'IndexSet( abs_p_interp_order, 5 )\n'); fprintf(fid,'IndexSet( abs_t_interp_order, 7 )\n'); fprintf(fid,'IndexSet( abs_nls_interp_order, 5 )}\n'); fclose(fid); end Q.INCLUDES= {fullfile(workfolder,'include.arts')}; ind=find(Q.SENSOR_LOS(:,2)>180); if length(ind)>0 Q.SENSOR_LOS(ind,2)=Q.SENSOR_LOS(ind,2)-360; end if Q.CLOUDBOX_DO==0 %use clearsky mode Q = asg2q( asgD, G, Q, workfolder ); Q.ABS_NLS=[]; Q = qarts_abstable( Q,9,100,[],[] ); [y1,ydata,dy] = arts_y( Q, workfolder ); else %do MC simulations y1=zeros(length(Q.F_GRID)*length(Q.SENSOR_LOS),1); std_err=unique(Q.CLOUDBOX.METHOD_PRMTRS.STD_ERR); %loop over different part of the "antenna" Q1=Q; for i=1:length(std_err) i ind=find(Q.CLOUDBOX.METHOD_PRMTRS.STD_ERR==std_err(i)); Q1.CLOUDBOX.METHOD_PRMTRS.STD_ERR=std_err(i); if i==1 Q1 = asg2q( asgD, G, Q1, workfolder ); Q1.ABS_NLS=[]; Q1 = qarts_abstable( Q1,9,100,[],[] ); end Q1.SENSOR_POS=Q.SENSOR_POS(ind,:); Q1.SENSOR_LOS=Q.SENSOR_LOS(ind,:); [y,ydata,dy] = arts_y( Q1, workfolder ); for j=1:length(ind) y1((ind(j)-1)*length(Q.F_GRID)+1:ind(j)*length(Q.F_GRID))=... y((j-1)*length(Q.F_GRID)+1:j*length(Q.F_GRID)); end end end end %convert to intensity boltzmann = constants('BOLTZMANN_CONST'); speed_light = constants('SPEED_OF_LIGHT'); f_grid=zeros(size(y1)); for j=1:length(Q.F_GRID) f_grid(j:length(Q.F_GRID):end)=Q.F_GRID(j); end %convert to Planck I=2*f_grid.^2*boltzmann/(speed_light^2).*y1; Y= i2planckTb(I,f_grid); %weight pencil beam simulations with antenna and sideband response Y=sensor.ch_resp*sensor.ant_resp*Y; return %----------------------------------------------------------------- function [Q,G]=asg_set_cb(Q,G) iwc_ind = find( strcmp( lower({G.NAME}), 'iwc field' ) ); [C1,C2,C3,C4,C5,C6]=deal([]); for i=1:size(G(iwc_ind).DATA,3) [c1,c2]=find(G(iwc_ind).DATA(:,:,i)); if ~isempty(c1) C1=min([C1,vec2col(c1)']); C2=max([C2,vec2col(c1)']); end if ~isempty(c2) C3=min([C3,vec2col(c2)']); C4=max([C4,vec2col(c2)']); end end for i=1:size(G(iwc_ind).DATA,1) [c1,c2]=find(squeeze(G(iwc_ind).DATA(i,:,:))); if ~isempty(c2) C5=min([C5,vec2col(c2)']); C6=max([C6,vec2col(c2)']); end end % if ~isempty(C1) %if C1 is empty we don't have any cloud particles p1 = G(iwc_ind).GRID1(C1); p2 = G(iwc_ind).GRID1(C2); lat1 = G(iwc_ind).GRID2(C3); lat2 = G(iwc_ind).GRID2(C4); z = p2z_cira86([p1 p2]',mean([lat1 lat2]),160); lon1 = G(iwc_ind).GRID3(C5); lon2 = G(iwc_ind).GRID3(C6); % Q.CLOUDBOX.LIMITS = [z(1)-2e3 z(2)+2e3 lat1 lat2 lon1 lon2]; Q.CLOUDBOX_DO=1; else %remove all particles in G ind =find(strncmp({G.NAME},'Particles',9)==0); G = G(ind); end return %---------------------------------------- function P=asg_atm_pb(G,Q,sensor) addpath /home/bengt/ARTS/atmlab/physics/thermodynamic p = asg_pathiwc( asgD, G, Q, 21e3 ); zvec=[0:250:18e3]; warning off for i=1:length(p) iwc(:,i)=interp1(p(i).z,p(i).iwc,zvec); lwc(:,i)=interp1(p(i).z,p(i).lwc,zvec); iwp(i)=p(i).iwp; rhi(:,i)=interp1(p(i).z,p(i).rhi,zvec); t(:,i)=interp1(p(i).z,p(i).t,zvec); end warning on %weight the individual profiles w=sensor.ant_resp(1,1:length(Q.F_GRID):end); P.z=zvec; P.iwc=iwc*w'; P.lwc=lwc*w'; P.iwp=iwp*w'; P.rhi=rhi*w'; P.t=t*w'; return function H = amsu_antenna(za_grid,aa_grid) Q = qarts; Q.SENSOR_DO = 1; Q.STOKES_DIM = 1; Q.ANTENNA_DIM = 2; Q.ATMOSPHERE_DIM = 3; Q.OUTPUT_FILE_FORMAT = 'ascii'; Q.F_GRID = ... fullfile( atmlab('ARTS_INCLUDES'), 'amsu', 'amsub.f_grid_fast.xml' ); Q.MBLOCK_ZA_GRID = za_grid; Q.MBLOCK_AA_GRID = aa_grid; Q.SENSOR_RESPONSE = qartsSensor; Q.SENSOR_RESPONSE.SENSOR_NORM = 1; Q.SENSOR_RESPONSE.ANTENNA_LOS = [0 0]; Q.SENSOR_RESPONSE.ANTENNA_RESPONSE = fullfile( atmlab('ARTS_INCLUDES'), ... 'amsu', 'amsub.antenna.xml' ); Q = arts_sensor( Q ); H = Q.SENSOR_RESPONSE; return %------- function [mc_std]=mc_error(w,std,std_max) mc_std=sqrt(std^2./w); mc_std(find(mc_std>std_max))=std_max; ind1=find(mc_std(:)==std_max); ind1b=ind1; ind2=find(mc_std(:)0 varp=std^2-sum(w(ind1b).^2.*mc_std(ind1b).^2); mc_std(ind2)=min(sqrt(varp/sum(w(ind2))./w(ind2)),std_max); ind1=find(mc_std(ind2)>std_max); if length(ind1)>0 ind1b=sort([ind1b,ind2(ind1)]) end ind2=setdiff(1:length(mc_std(:)),ind1b); end return