merge

README.md

@@ -0,0 +1,4 @@
+transcar
+========
+
+Blelly and Lilensten Transcar 1D ionosphere precipitation model

dir.transcar.server/BT_E1E2prev.csv

@@ -0,0 +1,33 @@
+52.726,63.557,49.521,59.731
+63.557,76.487,59.731,71.92
+76.487,91.921,71.92,86.47
+91.921,110.34,86.47,103.84
+110.34,132.34,103.84,124.57
+132.34,158.59,124.57,149.32
+158.59,189.93,149.32,178.86
+189.93,227.34,178.86,214.13
+227.34,272.,214.13,256.23
+272.,325.3,256.23,306.48
+325.3,388.94,306.48,366.46
+388.94,464.9,366.46,438.07
+464.9,555.57,438.07,523.55
+555.57,663.81,523.55,625.58
+663.81,793.02,625.58,747.38
+793.02,947.25,747.38,892.78
+947.25,1131.4,892.78,1066.3
+1131.4,1351.1,1066.3,1273.5
+1351.1,1613.5,1273.5,1520.8
+1613.5,1926.7,1520.8,1816.1
+1926.7,2300.5,1816.1,2168.5
+2300.5,2746.8,2168.5,2589.2
+2746.8,3279.5,2589.2,3091.3
+3279.5,3915.4,3091.3,3690.8
+3915.4,4674.4,3690.8,4406.3
+4674.4,5580.6,4406.3,5260.5
+5580.6,6662.2,5260.5,6280.2
+6662.2,7953.4,6280.2,7497.3
+7953.4,9494.7,7497.3,8950.3
+9494.7,11335.,8950.3,10685.
+11335.,13531.,10685.,12755.
+13531.,16152.,12755.,15227.
+16152.,19282.,15227.,18177.

dir.transcar.server/dir.input/DATCAR

@@ -1,16 +1,16 @@
 2					(kiappel, INTEGER)
 90kmmaxpt123.dat                                input file (initial ionospheric conditions)
 1.					time step (seconds) (REAL)
-30					number of seconds between fluid code outputs (seconds, REAL)
-2004022					date of simulation (YYYYDDD)
-21.					UT start time for simulation (seconds)
-1200.					duration of run (seconds)
+1.				number of seconds between fluid code outputs (seconds, REAL)
+2013090					date of simulation (YYYYDDD)
+32400.					UT start time for simulation (seconds)
+60.					duration of run (seconds)
 2					jpreci, precipation type (see table below, INTEGER)
-78.,16.    				location of simulation, geodetic latitude,longitude
+65.12,-147.43    				location of simulation, geodetic latitude,longitude
 0.					duration of convection (in seconds) before the reference (<= 0, we do not follow the field lines)  (tempsconv_1, REAL)
 0.					duration of convection (in seconds) before the reference (<= 0, we do not follow the field lines)  (tempsconv, REAL)
 0					time interval (seconds) between two tubes (step,REAL)
-30					time step between two calls to kinetic code (REAL)
+1.					time step between two calls to kinetic code (REAL)
 0.					transport (m / s) induced along the field line (vparaB, REAL)
 126.0					f10.7 index
 107.6					f10.7 3-month average centered on current month
@@ -25,8 +25,8 @@
 precinput.dat   	                electron precipitation distribution file (energy (eV) in first column, number flux (cm^-2 s^-1 eV^-1 sr^-1) in second column)
 1					precint, precipation energy interpolation (see table below, INTEGER)
 1					precext, precipation energy extrapolation (see table below, INTEGER)
-500.					precipitation start time (seconds)
-1000.					precipitation end time (seconds)
+32420.					precipitation start time (seconds from day start)
+32440.					precipitation end time (seconds from day start)
 
 
 

dir.transcar.server/dir.input/DATCAR~

@@ -0,0 +1,46 @@
+2					(kiappel, INTEGER)
+90kmmaxpt123.dat                                input file (initial ionospheric conditions)
+1.					time step (seconds) (REAL)
+30	.				number of seconds between fluid code outputs (seconds, REAL)
+2013090					date of simulation (YYYYDDD)
+32400.					UT start time for simulation (seconds)
+60.					duration of run (seconds)
+2					jpreci, precipation type (see table below, INTEGER)
+65.12,-147.43    				location of simulation, geodetic latitude,longitude
+0.					duration of convection (in seconds) before the reference (<= 0, we do not follow the field lines)  (tempsconv_1, REAL)
+0.					duration of convection (in seconds) before the reference (<= 0, we do not follow the field lines)  (tempsconv, REAL)
+0					time interval (seconds) between two tubes (step,REAL)
+1.					time step between two calls to kinetic code (REAL)
+0.					transport (m / s) induced along the field line (vparaB, REAL)
+126.0					f10.7 index
+107.6					f10.7 3-month average centered on current month
+63.6					ap index
+25.					convection electric field (mV/m)
+1.					O correction factor (REAL)
+1.					N2 correction factor (REAL)
+1.					O2 correction factor (REAL)
+1.					N correction factor (REAL)
+1.					H correction factor (REAL)
+-1.e-4					topside electron heat flux (mW/m^2, REAL)
+precinput.dat   	                electron precipitation distribution file (energy (eV) in first column, number flux (cm^-2 s^-1 eV^-1 sr^-1) in second column)
+1					precint, precipation energy interpolation (see table below, INTEGER)
+1					precext, precipation energy extrapolation (see table below, INTEGER)
+32420.					precipitation start time (seconds from day start)
+32440.					precipitation end time (seconds from day start)
+
+
+
+jpreci = 0 if sun only
+       = 1 if electron precipitation only
+       = 2 if sun + electron precipitation
+       = 3 if proton precipitation only
+       = 4 if proton and electron precipitation
+       = 5 if sun + proton precipitation
+       = 6 if proton, electron precipitation + sun
+
+precint = 0 log-linear interpolation
+	= 1 sample and hold interpolation
+
+precext = 0 log-linear extrapolation
+	= 1 truncation of distribution past energy range of input file
+	= 2 hold energy dist constant past energy range of input file

dir.transcar.server/loop_beams.sh

@@ -0,0 +1,32 @@
+#!/bin/bash
+# Michael Hirsch 2014
+# prototype for parallel execution, by looping 
+
+BeamLooper()
+{
+RTdir=$1
+[[ ! -d $RTdir ]] && { echo "Error, root directory $RTdir does not exist. Aborting"; exit 1; }
+
+Efile=$2;
+[[ ! -a $Efile ]] && { echo "Error, energy CSV file $Efile does not exist. Aborting"; exit 2; }
+
+
+cIFS=$IFS #we will restore this after the while loop
+IFS=,
+while read E1 E2 pr1 pr2 
+do
+	./run_beams.sh $RTdir $E1 $E2 $pr1 $pr2 # this is what parallel will do
+
+     #echo $RTdir $E1 $E2 $pr1 $pr2 #debug
+     #error handling
+     LastErr=$?
+     case $LastErr in 
+     99) echo "loop_beams: Debug exit $E1, code $LastErr"; exit 99 ;;
+     *) echo "loop_beams: $E1 returned $LastErr" ;;
+     esac
+
+done < $Efile
+IFS=$cIFS
+}
+
+BeamLooper $1 $2

dir.transcar.server/run_beams_parallel.sh

@@ -5,7 +5,11 @@
 # respective beam energies.
 # this program is meant to be run from dir.transcar.server directory- cd there first
 
+<<<<<<< HEAD
 BeamEnergyTableFN=Prosp8Beam
+=======
+BeamEnergyTableFN=BT_E1E2prev.csv
+>>>>>>> eb95003423a2db072bf674fa9aeb450e184846d3
 
 runBeams()
 {
@@ -16,11 +20,19 @@ E2=$2
 pr1=$3
 pr2=$4
 
+<<<<<<< HEAD
 RODIR=../../prosp0
 
 CurrDir="$RODIR/beam$E1"
 # freshen simulation directory
 [[ -d $RODIR ]] && { rm -rv $CurrDir; mkdir -v $CurrDir; }
+=======
+RODIR=../../iter0
+
+CurrDir="$RODIR/beam$E1"
+# freshen simulation directory
+[[ -d $RODIR ]] && { \rm -rv $CurrDir; mkdir -v $CurrDir; }
+>>>>>>> eb95003423a2db072bf674fa9aeb450e184846d3
 
 BMlog=$RODIR/Beams.log
 TCconfig=dir.input/DATCAR
@@ -51,10 +63,17 @@ PrecFN="$CurrDir/dir.input/precinput.dat"
   flux=$(echo "scale=4; $flux0 / 0.5 / $Esum / $dE" | bc) 
 
   #get 'padding' so that beams are truncated correctly in energy
+<<<<<<< HEAD
   pr=$pr1 #previous?
   dElow=$(echo "$E1 - $pr" | bc)
   Elow=$(echo "scale=4; $E1 - 0.5 * $dElow" | bc)
   ne=$pr2  #next?
+=======
+  pr=$pr1
+  dElow=$(echo "$E1 - $pr" | bc)
+  Elow=$(echo "scale=4; $E1 - 0.5 * $dElow" | bc)
+  ne=$pr2 
+>>>>>>> eb95003423a2db072bf674fa9aeb450e184846d3
   dEhigh=$(echo "$E2 - $ne" | bc)
   Ehigh=$(echo "scale=4; $E2 - 0.5 * $dEhigh" | bc)
 

optical_inversions/ISR_proc/barker_analyze.m

@@ -0,0 +1,23 @@
+isrz=100:10:310;
+isrneraw=[.5 1.4 2.4 3.8 4 3.9 3.8 3.5 3.3 3.4 3.7 3.8 3.8 3.85 3.85 3.5 3.25 2.6 2.4 2.3 2.2 2.1]*10^11;
+
+load barker.mat;            
+
+tind=min(find(tbark>=1+3/60));
+te_tiavg=mean(te_ti(:,tind:tind+5),2);            %use average Te/Ti to make up for the fact that the scanning took 2 minutes
+te_tibark=interpolate(te_tiavg,zbark,isrz,'lin','lin');
+nesimavg=mean(nesim(:,tind:tind+5),2);
+
+correction_fact=(te_tibark+1)/2;        %multiplicative correction factor
+nereal=correction_fact.*isrneraw;       %correct ne
+
+FS=16;
+plot([isrneraw(:),nereal(:)],isrz(:),'LineWidth',2);
+set(gca,'FontSize',FS)%,'FontWeight','bold');
+ax=axis;
+axis([ax(1:3),375])
+xlabel('n_e (m^{-3})');
+ylabel('altitude (km)');
+title('Barker-coded data from 19 Nov. 2001');
+legend('n_{e} with T_e/T_i=1','n_{e} using estimated T_e/T_i','Location','NorthWest');
+print -depsc barker_corr.eps

optical_inversions/ISR_proc/isdiff.m

@@ -0,0 +1,74 @@
+%CALCULATE CORRECTIONS TO BARKER-CODED DATA
+if ~exist('neplotbcorr')
+    isplot_barker_corr_wsim;
+end
+
+
+%SMOOTH THE BARKER PROFILES TO MAKE IT EASIER TO FIND MAXIMA
+mawindow=15;
+for k=1:length(tplotb)
+   neplotbsmooth(:,k)=smooth(neplotb(:,k),mawindow);
+   neplotbcorrsmooth(:,k)=smooth(neplotbcorr(:,k),mawindow);
+end
+for k=1:length(tplotf)
+   neplotfsmooth(:,k)=smooth(neplotf(:,k),4);
+end
+
+
+
+%FIND NMF2 AND HMF2 FOR EACH SIMULATION
+isearch=find(zplot>=175 & zplot<=400);
+for k=1:length(tplotb)
+   [nmf2b(k),imaxb]=max(neplotbsmooth(isearch,k));
+   hmf2b(k)=zplot(isearch(imaxb));
+end
+
+for k=1:length(tplotb)
+   [nmf2c(k),imaxc]=max(neplotbcorrsmooth(isearch,k));
+   hmf2c(k)=zplot(isearch(imaxc));
+end
+
+for k=1:length(tplotf)
+   [nmf2f(k),imaxf]=max(neplotfsmooth(isearch,k));
+   hmf2f(k)=zplot(isearch(imaxf));
+end
+
+
+
+%INTERPOLATE FITTED DATA ONTO BARKER TIME GRID
+nmf2fplot=interpolate(nmf2f',tplotf',tplotb','lin','lin')
+errb=abs(100*(10.^nmf2b-10.^nmf2fplot)./10.^nmf2fplot);
+errbcorr=abs(100*(10.^nmf2c-10.^nmf2fplot)./10.^nmf2fplot);
+itmean=min(find(tplotb>=1+2/60)):length(tplotb);
+avgerrb=mean(errb(itmean)); avgerrbcorr=mean(errbcorr(itmean));
+
+
+
+%PLOTS
+figure;
+
+subplot(211);
+plot(tplotb,10.^nmf2b,tplotf,10.^nmf2f,tplotb,10.^nmf2c,'LineWidth',2);
+axtemp=axis;
+ax=[min(tplotb),max(tplotb),2.5e11,8e11];
+axis(ax);
+set(gca,'FontSize',12);
+l=legend('raw n_e','fitted n_e','T_e/T_i adj. n_e','Location','NorthWest','Orientation','Horizontal');
+set(l,'FontSize',12);
+xlabel('UT (hrs)');
+ylabel('N_mF_2');
+set(gca,'XTick',xloc,'XTickLabel',xlab);
+
+subplot(212);
+plot(tplotb,errb,tplotb,errbcorr,'LineWidth',2);
+axtemp=axis;
+ax=[min(tplotb),max(tplotb),0,70];
+axis(ax);
+set(gca,'FontSize',12);
+l=legend(sprintf('raw n_e.  Avg:  %3.1f %%',avgerrb),sprintf('T_e/T_i adj. n_e.  Avg:  %3.1f %%',avgerrbcorr),'Location','NorthWest','Orientation','Horizontal');
+set(l,'FontSize',12);
+xlabel('UT (hrs)');
+ylabel('% Error in N_mF_2');
+set(gca,'XTick',xloc,'XTickLabel',xlab);
+
+print -depsc error_barker.eps

optical_inversions/ISR_proc/isplot_barker.m

@@ -0,0 +1,65 @@
+%GET THE PROCESSED ISR DATA
+load ./isrdata_barker.mat;
+
+
+
+%TICK MARKS FOR TIMES AXIS
+iminx=min(find(istimeb>=1+1/60));
+xloc=istimeb(iminx):2/60:istimeb(length(istimeb));
+for k=1:length(xloc)
+    toi=xloc(k);
+    labhr=int2str(floor(toi));
+    labmin=int2str(max(0,round(toi*60-str2num(labhr)*60)));
+    if length(labmin)<2
+        labmin=['0',labmin];
+    end
+    labsec=int2str(max(0,round(toi*3600-str2num(labmin)*60-str2num(labhr)*3600)));
+    if length(labsec)<2
+        labsec=['0',labsec];
+    end
+    xlab{k}=[labhr,':',labmin];
+end
+
+
+
+%GET RID OF SHITTY DATA POINTS
+ibad=find(isneb==0);
+igood=find(isneb(:)~=0);
+mindat=min(isneb(igood));
+isneb(ibad)=1e-6*mindat;
+
+
+
+%INTERPOLATE FOR IMAGE SCALING
+zinds=find(iszb>=0 & iszb<=600);
+minz=min(iszb(zinds)); maxz=max(iszb(zinds)); zspan=maxz-minz;
+
+tplot=istimeb;
+zplot=100:5:600;
+for k=1:length(tplot);
+   neplot(:,k)=interpolate(isneb(zinds,k),iszb(zinds),zplot,'loglin','loglin') ;
+end
+
+
+
+%PLOT AS A RASTERIZED IMAGE W/ SCALABLE FONTS
+figure;
+set(gcf,'PaperPosition',[0,0,8.5,3]);
+FS=12;
+h=imagesc(tplot,zplot,flipud(neplot),[10 11.8]);
+c=colorbar;
+h=gca;
+set(h,'FontSize',FS)
+set(h,'Xtick',xloc);
+set(h,'XtickLabel',xlab);
+yloc0=get(h,'YTick');
+dyloc=yloc0-round(minz);
+newyloc=round(maxz)-dyloc;
+set(h,'YTick',sort(newyloc));
+set(h,'YTickLabel',fliplr(yloc0));
+xlabel('UT (hrs)');
+ylabel('altitude (km)');
+%colormap('gray')
+ax=axis;
+ylabel(c,'log_{10}(raw n_{e}) (m^{-3})')
+print -depsc ne_barker_img.eps;