Cobweb plot

While conducting multi variable analysis, it may be useful to see how the output changes due to change in inputs. For example humidity, temperature and water-level may change the deflection of dam. Cobweb plot combined with normal distribution will help to visualize the risks easily. An example of such plot is shown below.

An octave (or matlab) code to generate such plot is given below. Basically, it requires to plot each entity separatly. The data should be normalized to set equal scale.

%shaded plot
clear all;clf; % clear figure (just in case something was there before)

epdata=[-0.0008959	-0.00059726	-0.0013775	-0.00080342	-0.002345	-0.0012454	-0.0014801	-0.00097711	-0.00096319	-0.00030882	-0.0015608	-0.00099932	-0.00082923	-0.00105	-0.0013915	-0.00021562	-0.0011739	-0.0015718	-0.00087249	-0.0016794	-0.00036352	-0.00059309	-0.0011483	-0.0010326	-0.00034531	-0.0020304	-0.0014554	-0.00092798	-0.0012557	-0.0016694	-0.0010306	-0.001245	-0.0013191	-0.0016508	-0.0010848	-0.0018287	-0.0010912	-0.0014441	-0.0012839	-0.00089405	-0.0024034	-0.0015758	-0.00026981	-0.0014717	-0.0023293	-0.0020914	-0.00076314	-0.00056843	-0.0014572	-0.00083516	-0.0012168	-0.0013884	-0.0012994	-0.0010078	-0.00045422	-0.0015285	-0.00061293	-0.00091066	-0.0014499	-0.0010975	-0.001501	-0.0012706	-0.0019573	-0.0019029	-0.0013711	-0.0019286	-0.00099036	-0.0016533	-0.001049	-0.0015944	-0.0016858	-0.0014946	0.000314	-0.0016982	-0.0012176	-0.0016377	-0.0012109	-0.0015727	-0.001435	-0.0010661	-0.0012172	-0.0020125]
tcdata=[31.007	17.279	71.559	42.942	30.757	64.937	56.682	34.429	48.118	51.053	46.382	70.517	49.346	56.116	38.708	39.931	61.857	52.765	58.071	55.669	34.416	38.345	36.623	42.745	66.059	83.693	40.504	62.218	64.917	20.47	65.233	37.162	69.894	27.246	3.5394	65.416	65.68	43.018	95.189	59.176	59.944	56.573	53.296	38.949	41.054	39.729	28.056	43.599	59.708	46.291	49.449	48.566	57.104	51.747	68.345	48.032	63.941	41.263	26.673	48.541	31.849	56.63	39.846	49.557	80.801	37.334	37.459	36.248	57.355	39.958	32.182	24.017	75.422	61.597	43.402	11.902	58.988	42.176	17.418	44.48	56.916	41.854]
tldata=[61.235	77.001	86.213	75.22	43.764	67.211	76.459	52.866	77.028	54.373	34.241	90.087	80.787	60.298	75.159	107.12	84.327	39.385	104.76	66.38	62.929	42.629	88.112	83.559	59.862	85.192	50.859	55.306	66.787	56.848	82.743	95.426	112.39	44.012	70.265	54.671	76.631	76.642	88.841	39.25	56.877	60.893	104.51	59.836	63.706	71.222	63.624	32.153	70.074	79.584	72.205	72.182	73.233	67.998	80.533	40.931	63.258	42.223	79.662	81.851	52.266	80.383	83.354	57.363	93.174	64.228	87.165	93.264	79.508	67.241	98.564	86.5	37.035	59.61	64.813	68.867	72.937	81.84	61.348	92.968	57.552	64.609]
dpdata=[-0.016616	0.0012712	-0.018897	-0.01016	-0.092141	-0.024361	-0.023951	-0.026684	-0.013089	-0.0056231	-0.065296	-0.011672	-0.0096546	-0.023527	-0.018468	0.0022514	-0.015527	-0.056125	-0.0052227	-0.0348	-0.0046483	-0.020469	-0.0072833	-0.010084	-0.0052599	-0.030438	-0.044018	-0.022032	-0.024795	-0.031457	-0.014002	-0.0055563	-0.0096733	-0.063971	0.003394	-0.046408	-0.017231	-0.019991	-0.016868	-0.028567	-0.061399	-0.036514	0.00094998	-0.034091	-0.050397	-0.03461	-0.011083	-0.023391	-0.027504	-0.0096521	-0.020398	-0.023604	-0.022343	-0.018593	-0.0046534	-0.055475	-0.011241	-0.03279	-0.0084095	-0.013231	-0.043265	-0.018161	-0.020591	-0.04901	-0.016887	-0.038989	-0.0064535	-0.0089957	-0.014825	-0.029465	-0.0045653	-0.0031735	0.013363	-0.040035	-0.024299	-0.00078227	-0.020854	-0.018003	-0.015243	-0.0074379	-0.02939	-0.042169]
% normalize
ep=epdata/min(epdata);
tc=tcdata/max(tcdata);
tl=tldata/max(tldata);
dp=dpdata/min(dpdata);
x=[0 1 2 3.5]*7;
x1=ones(1,length(ep))*x(1);
x2=ones(1,length(ep))*x(2);
x3=ones(1,length(ep))*x(3);
x4=ones(1,length(ep))*x(4);

%calculate mean
m_ep=mean(ep);
m_tc=mean(tc);
m_tl=mean(tl);
m_dp=mean(dp);
%standard deviation
sd_ep=std(ep);
sd_tc=std(tc);
sd_tl=std(tl);
sd_dp=std(dp);

%coords to plot normal distribution
nd_x_ep=[m_ep-4*sd_ep:sd_ep/20:m_ep+4*sd_ep]
nd_y_ep=1/sd_ep/(2*pi())^0.5*exp((nd_x_ep-m_ep).^2/-2/sd_ep^2)

nd_x_tc=[m_tc-4*sd_tc:sd_tc/20:m_tc+4*sd_tc]
nd_y_tc=1/sd_tc/(2*pi())^0.5*exp((nd_x_tc-m_tc).^2/-2/sd_tc^2)

nd_x_tl=[m_tl-4*sd_tl:sd_tl/20:m_tl+4*sd_tl]
nd_y_tl=1/sd_tl/(2*pi())^0.5*exp((nd_x_tl-m_tl).^2/-2/sd_tl^2)

nd_x_dp=[m_dp-4*sd_dp:sd_dp/20:m_dp+4*sd_dp]
nd_y_dp=1/sd_dp/(2*pi())^0.5*exp((nd_x_dp-m_dp).^2/-2/sd_dp^2)

figure 1
% plot vertical lines
hold on;
plot (x1,ep,"o"); 
plot (x2,tc,"o"); 
plot (x3,tl,"o"); 
plot (x4,dp,"o"); 
%plot connecting lines
p=[ep; tc; tl; dp ];
plot (x, p, "color",[.8 .8 .8] );

%plot the normal disribution
patch(nd_y_ep*2+x(1),nd_x_ep,'facecolor',[1 1 0.8],'edgecolor',[1 0 0], 'facealpha',0.3)
patch(nd_y_tc*2+x(2),nd_x_tc,'facecolor',[1 1 0.8],'edgecolor',[1 0 0], 'facealpha',0.3)
patch(nd_y_tl*2+x(3),nd_x_tl,'facecolor',[1 1 0.8],'edgecolor',[1 0 0], 'facealpha',0.3)
patch(-nd_y_dp*2+x(4),nd_x_dp,'facecolor',[1 1 0.8],'edgecolor',[0 0 1], 'facealpha',0.3,"linewidth", 3)

%yaxis lines
ymin=-0.5;ymax=1.5;
plot ([x(1) x(1)],[ymin ymax],"color",[0 0 0] );
plot ([x(2) x(2)],[ymin ymax],"color",[0 0 0] );
plot ([x(3) x(3)],[ymin ymax],"color",[0 0 0] );
plot ([x(4) x(4)],[ymin ymax],"color",[0 0 0] );
%other properties
ylim([ymin ymax])
axis off;
set(gca, "linewidth", 2, "fontsize", 12);
%labels
%xlabels
text(x(1)+0.3,ymin,"Strain","color","black","fontsize", 18)
text(x(2)+0.3,ymin,"TC","color","black","fontsize", 18)
text(x(3)+0.3,ymin,"TL","color","black","fontsize", 18)
text(x(4)-3.5,ymin,"Max. Deflection","color","black","fontsize", 18)
%informations
text(x(1)+0.3,ymax-.2,["Mean=" num2str(-0.0012, "%5.2e") char(10) "SD=" num2str(0.005, "%5.2e")],"color","black","fontsize", 15)
text(x(2)+0.3,ymax-.2,["Mean=" num2str(50, "%5.0f") char(10) "SD=" num2str(15, "%5.0f") " years"],"color","black","fontsize", 15)
text(x(3)+0.3,ymax-.2,["Mean=" num2str(70, "%5.0f") char(10) "SD=" num2str(20, "%5.0f") " years"],"color","black","fontsize", 15)
text(x(4)-3.5,ymax-.2,["Mean=" num2str(mean(dpdata), "%5.2e") char(10) "SD=" num2str(std(dpdata), "%5.2e") " m"],"color","black","fontsize", 15)

hold off;
print -dpng figure1.png

Evolution and entropy

Some salt required.

Once born we slowly evolve physically by gaining tissues and bones. The food we eat is transformed into the machinaries of body. This indicates that the entropy or randomness decreases over time to create a useful thing. But once we cross certain age decay starts and the biological machinaries starts to disfunction. It seems decay and dissolution is the final destiny. And in long term degeneration is imenent.

If so how about the evolution theory? Why should a simple being evolve to a complex one? More natural way would be a complex being degenerate into simple ones. Because this is the way we observe, its easy to break things than to create one. May be we human too are degenerated form of more complex superhumans, not the complex form of monkeys.

Code aster tutorial for load varying with time and depth

Here I post two tutorial in code aster/ salome meca as follows:
1. Applying varying load
In this tutorial, simple vertical element is considered in which a vertically varying load is applied by defining a function.

2. Analysis of dam for varying water level
This is a typical analysis of dam where water level is assumed to fluctuate with time. Since two variables changes i.e. time and water depth, we need to make an external file to read the pressure (load). This tutorial will show how to do it.

Google script to extract email information

This code helps to get summary of email such as subject, received date etc. Code is for the unread emails; remove lines with if condition check to read all emails.
Make a new google sheet> tools>script editor and copy and paste it. You will require to authorize the use of code while executing it.

function retriveUnreadMessages() {
// imports subject, from (email), received date of all unread email in inbox of gmail account
// the data is then written in the google spreadsheet
//reference:https://developers.google.com/apps-script/reference/spreadsheet/spreadsheet#appendRow(Object)
//

//-----------google sheet part
var ss=SpreadsheetApp.getActiveSheet();
var date1=new Date().toISOString().slice(0, 10);
ss.appendRow(["***** new execution = ",date1,"*****"]);// append a new information line where the new execution is done because data is appended at the end

//-----------------------------------------gmail part-----
// Log the subject lines of your Inbox
var threads = GmailApp.getInboxThreads();
// number of message to retrive
var n= threads.length;
for (var i = 0; i < n; i++) {

// check if there is unread messages in a thread
if(threads[i].isUnread)
{

var subject1=threads[i].getFirstMessageSubject(); // subjet of thread
var email1=threads[i].getMessages(); // messages inside the thread
// check if the individual message is unread
for (var j = 0; j < email1.length; j++) {
if(email1[j].isUnread)
{
var message_sender=email1[j].getFrom();
var message_sendDateTime=email1[j].getDate();
}
ss.appendRow([i,j,subject1,message_sender,message_sendDateTime]); //write data to spreadsheet.. it is appended to new row

}
}
}
}

Fractals in matlab/octave

Codes in matlab/octave for generating some popular fractals:

Mandelbrot

mandelbrot

function mandelbrotmain()
clear variables; clc;

maxit=200;
x=[-2,1];
y=[-1 1];

xpix=601;
ypix=401;

x=linspace(x(1),x(2),xpix);
y=linspace(y(1),y(2),ypix);
[xG, yG]=meshgrid(x,y);

c=mb(maxit,xG,yG);

figure
imagesc(x,y,c);
colormap([1 1 1;0 0 0]);
axis on;
grind on;

endfunction

function count=mb(maxItr,xG,yG)
c=xG+1i*yG;
count=ones(size(c));
z=c;

for n=1:maxItr
z=z.*z+c;
inside=abs(z)<=2;
count=count+inside;
endfor
endfunction

 

Sierpinski triangle
sierpinski

clc;
clear variables;
close windows;

clf

N=500;
x=zeros(1,N);y=x;r=x;
for a=2:N
c=randi([0 2]);
r(a)=c;
switch c
case 0
x(a)=0.5*x(a-1);
y(a)=0.5*y(a-1);
case 1
x(a)=0.5*x(a-1)+.25;
y(a)=0.5*y(a-1)+sqrt(3)/4;
case 2
x(a)=0.5*x(a-1)+.5;
y(a)=0.5*y(a-1);
end
end
plot(x,y,’o’)
title(‘Sierpinski’s triangle’)
legend(sprintf(‘N=%d Iterations’,N))

 

Fractal tree
fractaltrree

function treemain
tic
clear all;
depth = 7;
figure 1;
hold on;

drawTree2(0, 0, 90, depth);
toc;
endfunction

function drawTree2(x1, y1, angle, depth)
deg_to_rad = pi / 180.0;
rot=30; #degrees :: rotation from second iteration
branchLength=100*depth;
if (depth != 0)
x2 = x1 + cos(angle * deg_to_rad) * branchLength;
y2 = y1 + sin(angle * deg_to_rad) * branchLength ;
line([x1, x2], [y1, y2],’LineWidth’,depth);
drawTree2(x2, y2, angle – rot, depth – 1);
drawTree2(x2, y2, angle + rot, depth – 1);

endif
endfunction

 

A note on fractal (from HH Hardey):
There are upper and lower limits to describe natural objects by fractals. A figure may look like a fern leaf, but continued generation of the “leaf on larger and larger scales produces only a larger and larger fern leaf. A fern “plant” will never be produced. Similarly, there is a smallest scale to which a fractal description applies. As the overall fern leaf pattern is repeated on smaller and smaller scales, eventually the scale of a single fern plant cell will be reached. A single cell does not look like the overall shape.