Initial version from Martin

atovs_read_find.m

+%% ATOVS_READ_FIND   Read header and data entries of MHS data file. Find intrusion of Moon in DSV.
+%
+% The function extracts header information and selected data entries of an MHS data file. Then it searches for an intrusion of the Moon.
+%
+% FORMAT   [hdrinfo, data, err] = atovs_read_headerMHS19_all2( file_name );
+%
+% IN    file_name   Name of an MHS data file.
+%
+
+function [hdrinfo, data, err] = atovs_read_headerMHS19_all2(filename)
+mean_scannumber = 1000;
+ang_diam = 0.5;
+% we want to read this many bytes (MHS Data Set Header Record size)
+rec_len = 3072;
+
+% UNCOMPRESS if needed
+tmpdir = create_tmpfolder;
+c = onCleanup(@() rmdir(tmpdir,'s'));
+filename = uncompress(filename,tmpdir);
+if isempty(filename), error(errId,'Uncompressing failed'); end
+
+% open a file
+file_id = fopen( filename, 'r' );
+fseek( file_id, 0, -1 );
+firstbytes=fread( file_id, 3, 'uint8=>char' );
+if strcmp(firstbytes', 'NSS')
+    fseek( file_id, 0, -1 );
+else
+    fseek( file_id, 512, -1 );    
+end
+
+% file names not starting with NSS have 512 additional lines.
+
+% read header
+[header, count] = fread( file_id, rec_len, 'uint8' );
+err = header;
+if count ~= rec_len
+
+% close a file
+    fclose( file_id );    
+    disp( 'Error. Input file is not valid.' );
+    err = 1;
+    return
+end
+
+% Extracting header fields.
+% Index numbers correspond to start/end octet from the KLM User Guide.
+% http://webapp1.dlib.indiana.edu/virtual_disk_library/index.cgi/2790181/FID3711/klm/html/c8/s8317-2.htm
+hdrinfo.creationsite  = char(header(1:3)');
+hdrinfo.datasetname   = char(header(23:64)');
+hdrinfo.noaascid      = extract_uint16(header, 73, 74);
+hdrinfo.startyear     = extract_uint16(header, 85, 86);
+hdrinfo.startday      = extract_uint16(header, 87, 88);
+hdrinfo.nlines        = extract_uint16(header, 135, 136);
+
+% init data format
+[rec_format, rec_len, nchan, nfovs] = atovs_define_amsubl1c;
+
+% skip the header
+fseek( file_id, 3072, -1 );
+if ~strcmp(firstbytes', 'NSS') 
+    fseek( file_id, 512, 0 );
+end
+ 
+%file names not starting with NSS have 512 additional lines.
+rec_len = 3072;
+% read all records
+[record, count] = fread( file_id, uint32(rec_len) * uint32(hdrinfo.nlines), 'uint8' );
+% close a file
+fclose( file_id );
+% number of scan lines read
+nlines_read = count / rec_len;
+
+% if amount of data read is less than asked
+if count < rec_len * hdrinfo.nlines
+% if some scan lines are missing, we still can go on
+    if iswhole( nlines_read )
+        disp( 'Warning. Some scanlines are missing.' );
+% if number of scan lines is not integer, part of a record is missing
+    else
+        disp( 'Error. Input file is corrupt.' );
+        data = [];
+        err = 1;
+        return
+    end
+end
+
+for h=1:nlines_read
+    data.year(h)     = extract_uint16(record, 3, 4);
+    data.day(h)      = extract_uint16(record, 5, 6);
+    data.time(h)     = double(extract_uint32(record, 3072*(h-1)+9, 3072*(h-1)+12)/1000);
+    data.height(h)   = extract_uint16(record, 3072*(h-1)+211, 3072*(h-1)+212);
+    data.latfov45(h) = extract_int32(record, 3072*(h-1)+753+8*44, 3072*(h-1)+753+3+8*44);
+    data.lonfov45(h) = extract_int32(record, 3072*(h-1)+753+4+8*44, 3072*(h-1)+753+7+8*44);
+    data.latfov46(h) = extract_int32(record, 3072*(h-1)+753+8+8*44, 3072*(h-1)+753+11+8*44);
+    data.lonfov46(h) = extract_int32(record, 3072*(h-1)+753+12+8*44, 3072*(h-1)+753+15+8*44);
+    data.dsvcounts1(h) = extract_uint16(record, 3072*(h-1)+2569, 3072*(h-1)+2569+1);
+    data.dsvcounts2(h) = extract_uint16(record, 3072*(h-1)+2569+12, 3072*(h-1)+2569+13);
+    data.dsvcounts3(h) = extract_uint16(record, 3072*(h-1)+2569+24, 3072*(h-1)+2569+25);
+    data.dsvcounts4(h) = extract_uint16(record, 3072*(h-1)+2569+36, 3072*(h-1)+2569+37);
+    data.dsvcounts11(h) = extract_uint16(record, 3072*(h-1)+2569+2, 3072*(h-1)+2569+3); 
+    data.dsvcounts12(h) = extract_uint16(record, 3072*(h-1)+2569+4, 3072*(h-1)+2569+5);
+    data.dsvcounts13(h) = extract_uint16(record, 3072*(h-1)+2569+6, 3072*(h-1)+2569+7);
+    data.dsvcounts14(h) = extract_uint16(record, 3072*(h-1)+2569+8, 3072*(h-1)+2569+9);
+    data.dsvcounts15(h) = extract_uint16(record, 3072*(h-1)+2569+10, 3072*(h-1)+2569+11);
+    data.dsvcounts21(h) = extract_uint16(record, 3072*(h-1)+2569+14, 3072*(h-1)+2569+15); 
+    data.dsvcounts22(h) = extract_uint16(record, 3072*(h-1)+2569+16, 3072*(h-1)+2569+17);
+    data.dsvcounts23(h) = extract_uint16(record, 3072*(h-1)+2569+18, 3072*(h-1)+2569+19);
+    data.dsvcounts24(h) = extract_uint16(record, 3072*(h-1)+2569+20, 3072*(h-1)+2569+21);
+    data.dsvcounts25(h) = extract_uint16(record, 3072*(h-1)+2569+22, 3072*(h-1)+2569+23);
+    data.dsvcounts31(h) = extract_uint16(record, 3072*(h-1)+2569+26, 3072*(h-1)+2569+27); 
+    data.dsvcounts32(h) = extract_uint16(record, 3072*(h-1)+2569+28, 3072*(h-1)+2569+29);
+    data.dsvcounts33(h) = extract_uint16(record, 3072*(h-1)+2569+30, 3072*(h-1)+2569+31);
+    data.dsvcounts34(h) = extract_uint16(record, 3072*(h-1)+2569+32, 3072*(h-1)+2569+33);
+    data.dsvcounts35(h) = extract_uint16(record, 3072*(h-1)+2569+34, 3072*(h-1)+2569+35);
+    data.dsvcounts41(h) = extract_uint16(record, 3072*(h-1)+2569+38, 3072*(h-1)+2569+39); 
+    data.dsvcounts42(h) = extract_uint16(record, 3072*(h-1)+2569+40, 3072*(h-1)+2569+41);
+    data.dsvcounts43(h) = extract_uint16(record, 3072*(h-1)+2569+42, 3072*(h-1)+2569+43);
+    data.dsvcounts44(h) = extract_uint16(record, 3072*(h-1)+2569+44, 3072*(h-1)+2569+45);
+    data.dsvcounts45(h) = extract_uint16(record, 3072*(h-1)+2569+46, 3072*(h-1)+2569+47);
+end
+ 
+% Berechne Signal von jeder Deep-Space-View zum Zeitpunkt des Monddurchgangs in Kanal 1
+% Das gesuchte Maximum muss mindestens 120 Scans von Anfang und Ende des Erdumlaufs entfernt sein.
+intr_start = 121;
+intr_end = nlines_read-131;
+y1 = double(squeeze(data.dsvcounts11));
+y2 = double(squeeze(data.dsvcounts21));
+y3 = double(squeeze(data.dsvcounts31));
+y4 = double(squeeze(data.dsvcounts41));
+for i = intr_start:intr_end
+	A(1,i) = median(y1(i:i+10))-median(y1(i+60:i+70));
+	A(2,i) = median(y2(i:i+10))-median(y2(i+60:i+70));
+	A(3,i) = median(y3(i:i+10))-median(y3(i+60:i+70));
+	A(4,i) = median(y4(i:i+10))-median(y4(i+60:i+70));
+end
+M(1) = max(A(1,:));
+M(2) = max(A(2,:));
+M(3) = max(A(3,:));
+M(4) = max(A(4,:));
+maximum = max(max(M));
+[x,y]=find(A==maximum);
+longitude=(double(data.lonfov45(y))+double(data.lonfov46(y)))/2E4; % longitude of nadir at maximum of light curve
+latitude=(double(data.latfov45(y))+double(data.latfov46(y)))/2E4; % latitude of nadir at maximum of light curve
+height=double(data.height(y))/10; % altitude of satellite at maximum of light curve
+Z0 = ['Map Ref ',num2str(round(longitude,1)),'deg ',num2str(round(latitude,1)),'deg'];
+disp(Z0)
+Z1 = ['Altitude of Satellite = ',num2str(height),' km'];
+disp(Z1)
+Z2 = ['Maximum counts found in DSV no. ',num2str(x),' for scan ID ',num2str(y)];
+disp(Z2)
+Z3 = ['Date (yyyy-DDD   HR:MN) ',num2str(data.year(y)),'-',num2str(data.day(y)),'   ',num2str(fix(data.time(y)/3600)),':',num2str(fix(data.time(y)/60-60*fix(data.time(y)/3600)))];
+disp(Z3)
+end
+%big endian
+ function value = extract_uint16(bytearr, startbyte, endbyte);
+ value = typecast(uint8(bytearr(endbyte:-1:startbyte)), 'uint16');
+ end
+ 
+ %function value = extract_int16(bytearr, startbyte, endbyte);
+ %value = typecast(uint8(bytearr(endbyte:-1:startbyte)), 'int16');
+ %end
+ 
+ function value = extract_int32(bytearr, startbyte, endbyte);
+ value = typecast(uint8(bytearr(endbyte:-1:startbyte)), 'int32');
+ end
+ 
+ %function value = extract_uint32(bytearr, startbyte, endbyte);
+ %value = typecast(uint8(bytearr(endbyte:-1:startbyte)), 'uint32');
+ %end
+
+%Little endian
+%function value = extract_uint16(bytearr, startbyte, endbyte);
+%value = typecast(uint8(bytearr(startbyte:endbyte)), 'uint16');
+%end
+
+function value = extract_int16(bytearr, startbyte, endbyte);
+value = typecast(uint8(bytearr(endbyte:-1:startbyte)), 'int16');
+end
+
+%function value = extract_int32(bytearr, startbyte, endbyte);
+%value = typecast(uint8(bytearr(startbyte:endbyte)), 'int32');
+%end
+
+function value = extract_uint32(bytearr, startbyte, endbyte);
+value = typecast(uint8(bytearr(endbyte:-1:startbyte)), 'uint32');
+end