#!/usr/bin/perl -w
#
# makelattice.pl [a] [c] [rows] [columns]
#
#
use strict;

my @global_down = ();
my @global_up = ();
my @type = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1);
my $a = shift or die "Missing lateral lattice constant";
my $c = shift or die "Missing normal lattice constant";
my $m = shift or die "Missing number of rows";
my $n = shift or die "Missing number of columns";

init_lipids();


my $rho = 20.;
my $kn = 100.;
my $chin = 30.;
my $N = 16;

my $v_aa = -2. * ($kn + 3.) / $rho;
my $v_bb = 0.1;
my $v_ab = $chin / $rho + .5 * ($v_aa + $v_bb);
my $w_aaa = 1.5 * ($kn + 2.) / ($rho * $rho);
my $w_aab = $w_aaa;
my $w_abb = $w_aaa;
my $w_bbb = 0.0;
my $header = sprintf('v=%lg %lg %lg w=%lg %lg %lg %lg',
$v_aa, $v_ab, $v_bb, $w_aaa, $w_aab, $w_abb, $w_bbb);

my $lx = 50.;
my $ly = ($n) * $a;
my $lz = ($m) * $a * sqrt(3.); # nicht 100% sicher

print "# L=$lx $ly $lz n=". 4.*$m*$n ."  N=$N t=0\n";
print "# $header\n";

for (my $i = 0; $i < $n; $i++) {

	for (my $j = 0; $j < $m; $j++) {

		my $y0 = $i * $a;
		my $z0 = sqrt(3.) * $j * $a;

		my $y1 = $y0 + .5 * $a;
		my $z1 = $z0 + sqrt(.75) * $a;

		put_lipid_down(.5 * $lx - .5 * $c, $y0, $z0);
		put_lipid_down(.5 * $lx - .5 * $c, $y1, $z1);
		put_lipid_up(.5 * $lx + .5 * $c, $y0, $z0);
		put_lipid_up(.5 * $lx + .5 * $c, $y1, $z1);
	}
}

sub is_float { $_[0] =~ /^[+-]?\d+\.?\d*$/ } 

sub box_muller
{
	my $m = shift;
	my $s = shift;
	my $u1 = rand();
	my $u2 = rand();
	my $z;

#do {
		$z = sqrt(-2. * log(1. - $u1)) * cos(2. * 3.141592654 * $u2);
#	} while(!is_float($z));

	return $m + $z * $s; 
}


sub put_lipid_down
{
	my $x = shift;
	my $y = shift;
	my $z = shift;

	my $j = int(@global_up * rand());
	
	for (my $i = 0; $i < $N; $i++) {
		my $vx = box_muller(0, 1.);
		my $vy = box_muller(0, 1.);
		my $vz = box_muller(0, 1.);

		my $mx = $x + $global_down[$j][$i][0];
		my $my = $y + $global_down[$j][$i][1];
		my $mz = $z + $global_down[$j][$i][2];
		
		print "$mx $my $mz $vx $vy $vz $type[$i]\n";
	}
}

sub put_lipid_up
{
	my $x = shift;
	my $y = shift;
	my $z = shift;

	my $j = int(@global_up * rand());
	
	for (my $i = 0; $i < $N; $i++) {
		my $vx = box_muller(0, 1.);
		my $vy = box_muller(0, 1.);
		my $vz = box_muller(0, 1.);

		my $mx = $x + $global_up[$j][$i][0];
		my $my = $y + $global_up[$j][$i][1];
		my $mz = $z + $global_up[$j][$i][2];
		
		print "$mx $my $mz $vx $vy $vz $type[$i]\n";
	}
}

sub init_lipids
{
my @down = ();
my @up = ();
$down[0] = <<EOF;
24.4742 45.6689 27.3063
24.3322 45.6338 27.3406
24.0908 45.5269 27.289
23.8503 45.5688 27.2602
23.561 45.6763 27.2658
23.3192 45.6968 27.2928
23.0227 45.6027 27.1753
22.8319 45.6443 27.152
22.3567 45.613 27.2463
22.1138 45.6469 27.2582
21.8882 45.6915 27.1952
21.5551 45.773 27.2288
21.0883 45.7003 27.4629
20.8201 45.6087 27.7277
20.3294 45.5529 27.8429
19.9918 45.3707 28.1133
EOF

$down[1] = <<EOF;
25.3028 55.4412 35.7757
25.014 55.3242 35.8584
24.7135 55.3879 35.9033
24.5782 55.3587 35.8963
24.1739 55.2421 35.95660
23.7807 55.2501 35.9985
23.5914 55.226 35.955
23.4353 55.163 35.8983
23.0452 55.1738 35.9404
22.8602 55.1292 35.929
22.5032 55.0827 35.9554
22.2056 55.005 36.2134
21.5523 54.877 36.3874
21.1578 54.7505 36.808
21.0901 54.8272 36.94481
20.966 55.0634 37.0391
EOF

$down[2] = <<EOF;
24.9788 66.3087 32.139
24.7203 66.2508 32.02755
24.5162 66.4504 32.1166
24.4503 66.6933 32.253
24.2742 67.0688 32.2618
24.0552 67.225 32.4282
23.8382 67.2299 32.4811
23.4362 67.2501 32.57670
23.2412 67.1918 32.5477
23.0959 67.1447 32.5631
22.7679 67.2757 32.5317
22.4199 67.4684 32.3775
22.0294 67.802 32.2851
21.6965 68.2921 32.3924
21.665 68.4275 32.5469
21.7322 68.6555 32.7854
EOF

$down[3] = <<EOF;
24.9585 7.73382 1.31989
24.6706 7.74724 1.41321
24.37 7.72314 1.45795
24.1556 7.69011 1.42307
23.5834 7.80469 1.37649
23.3593 7.84233 1.44047
23.1566 7.91952 1.38301
22.8023 8.00979 1.56393
22.6103 7.97676 1.41939
22.2687 7.92795 1.41192
22.0004 7.90735 1.4554
21.6457 7.77213 1.63578
21.14 7.79814 1.7569
20.7158 7.80934 1.91439
19.9377 7.75307 1.63023
19.6344 7.71739 1.35361
EOF

$down[4] = <<EOF;
23.8317 9.73177 39.2379
23.818 9.87981 39.258
23.7649 10.1271 39.4394
23.748 10.1998 39.7462
23.8389 10.0634 40.22240
23.9248 10.2178 40.41710
24.16 10.3809 40.4019
24.2528 10.5352 40.4227
24.219 10.9441 40.6083
23.8203 11.0822 40.71651
23.5756 11.088 40.7071
23.1393 11.1591 40.8005
22.6115 11.203 40.7433
22.1609 10.9622 40.67721
22.0167 10.9675 40.5894
21.8326 10.6856 40.6276
EOF

$up[0] = <<EOF;
24.2815 37.3241 61.3428
24.71 37.111 61.6293
25.0996 37.0886 61.6921
25.5604 36.9767 61.7147
25.7598 36.9402 61.6467
26.1502 36.9501 61.672
26.5473 36.9483 61.6671
26.6372 36.9194 61.6827
26.9794 36.8423 61.7053
27.4354 36.9201 61.6886
27.5599 36.919 61.7156
27.95 36.8554 61.7315
28.2555 37.0078 62.1349
28.5567 37.4029 62.3879
28.6927 37.6134 62.425
28.8053 37.7548 62.3872
EOF

$up[1] = <<EOF;
24.8243 10.1177 3.08729
25.0527 10.0301 3.0354
25.4556 10.1105 3.01968
25.7038 10.1589 3.08565
26.0199 10.1225 3.11201
26.2084 10.1989 3.06321
26.4557 10.2003 3.05221
26.6679 10.0881 3.13275
26.889 10.0996 3.17461
27.4551 10.0819 3.16477
27.6144 10.1237 3.21539
27.8549 10.0982 3.15393
28.3374 10.0849 3.15592
28.7391 10.1936 3.27896
28.9302 10.1705 3.40599
29.1614 10.0623 3.78122
EOF

$up[2] = <<EOF;
24.9808 38.1899 56.3664
25.2515 38.2244 56.4042
25.5127 38.218 56.3554
25.8745 38.1792 56.4273
26.0823 38.2247 56.4333
26.4736 38.0384 56.5156
26.8675 38.1115 56.5852
27.0566 38.1615 56.5628
27.3527 38.1236 56.3104
27.4842 38.1129 56.2829
27.8351 38.0675 56.3426
28.1921 37.9984 56.4289
28.676 37.9296 56.4397
28.7536 37.8694 56.5083
29.3216 37.7704 56.8744
29.6199 37.6971 57.0518
EOF

	for (my $j = 0; $j < @down; $j++) {
		my @monomers = split /\n/, $down[$j];
		my @com;
		for (my $i = 0; $i < @monomers; $i++) {
			my @coord = split / /, $monomers[$i];
			$global_down[$j][$i][0] = $coord[0];
			$global_down[$j][$i][1] = $coord[1];
			$global_down[$j][$i][2] = $coord[2];
			$com[0] += $coord[0];
			$com[1] += $coord[1];
			$com[2] += $coord[2];
		}
		$com[0] /= @monomers;
		$com[1] /= @monomers;
		$com[2] /= @monomers;
		for (my $i = 0; $i < @monomers; $i++) {
			$global_down[$j][$i][0] -= $com[0];
			$global_down[$j][$i][1] -= $com[1];
			$global_down[$j][$i][2] -= $com[2];
		}
	}

	for (my $j = 0; $j < @up; $j++) {
		my @monomers = split /\n/, $up[$j];
		my @com;
		for (my $i = 0; $i < @monomers; $i++) {
			my @coord = split / /, $monomers[$i];
			$global_up[$j][$i][0] = $coord[0];
			$global_up[$j][$i][1] = $coord[1];
			$global_up[$j][$i][2] = $coord[2];
			$com[0] += $coord[0];
			$com[1] += $coord[1];
			$com[2] += $coord[2];
		}
		$com[0] /= @monomers;
		$com[1] /= @monomers;
		$com[2] /= @monomers;
		for (my $i = 0; $i < @monomers; $i++) {
			$global_up[$j][$i][0] -= $com[0];
			$global_up[$j][$i][1] -= $com[1];
			$global_up[$j][$i][2] -= $com[2];
		}
	}
}