text:
Determine Collision Point and Normal (3D SAT AABB-AABB)

Recommended Posts

AABB.cpp

int AABB::supportFaceCount()
{
// there are only three directions for every face of an AABB box.
return 3;
}

XMVECTOR AABB::supportFaceDirection(int i)
{
// the three axes of an AABB box. along the x, y and z axis.
static const XMVECTOR s_aabbAxes[] = { XMVectorSet(1, 0, 0, 0), XMVectorSet(0, 1, 0, 0), XMVectorSet(0, 0, 1, 0) };
return s_aabbAxes[i];
}

int AABB::supportEdgeCount()
{
// there are only three directions for every edges of an AABB box.
return 3;
}

XMVECTOR AABB::supportEdgeDirection(int i)
{
// every edge go along the x y, or z axis.
static const XMVECTOR s_aabbEdges[] = { XMVectorSet(1, 0, 0, 0), XMVectorSet(0, 1, 0, 0), XMVectorSet(0, 0, 1, 0) };
return s_aabbEdges[i];
}

void AABB::supportInterval(XMVECTOR direction, float& min, float& max)
{
XMVECTOR centre = XMVectorSet(Center[0], Center[1], Center[2], 1);
// projection of the box centre
float p = XMVector3Dot(centre, direction).m128_f32[0];

// projection of the box extents
float rx = fabs(direction.m128_f32[0]) * Radius[0];
float ry = fabs(direction.m128_f32[1]) * Radius[1];
float rz = fabs(direction.m128_f32[2]) * Radius[2];

// the projection interval along the direction.
float rb = rx + ry + rz;
min = p - rb;
max = p + rb;
}
bool ObjectsSeparatedAlongDirection(XMVECTOR& direction, AABB* a, AABB* b)
{
float mina, maxa;
float minb, maxb;
a->supportInterval(direction, mina, maxa);
b->supportInterval(direction, minb, maxb);
return (mina > maxb || minb > maxa);
}

bool ObjectsIntersected(AABB* a, AABB* b)
{
// test faces of A
for(int i = 0; i < a->supportFaceCount(); i++)
{
XMVECTOR direction = a->supportFaceDirection(i);
if(ObjectsSeparatedAlongDirection(direction, a, b))
return false;
}

// test faces of B
for(int i = 0; i < b->supportFaceCount(); i++)
{
XMVECTOR direction = b->supportFaceDirection(i);
if(ObjectsSeparatedAlongDirection(direction, a, b))
return false;
}

// test cross product of edges of A against edges of B.
for(int i = 0; i < a->supportEdgeCount(); i++)
{
XMVECTOR edge_a = a->supportEdgeDirection(i);

for(int j = 0; j < b->supportEdgeCount(); j++)
{
XMVECTOR edge_b = b->supportEdgeDirection(j);

XMVECTOR direction = XMVector3Cross(edge_a, edge_b);

if(ObjectsSeparatedAlongDirection(direction, a, b))
return false;
}
}
return true;
}

Share on other sites

You don't need to test edge directions for two AABB. This is only necessary for OBB.

For two AABB everything boils down to determine the two faces that are in contact. Then you can clip one against the side planes of the other. This gives you the contact points. For the normal you can choose either face direction.

I gave a talk about contact generation some time ago. It doesn't handle AABB vs AABB, but the more general OBB vs OBB. Your case is simpler since you don't need to handle the edge cases.

Edited by Dirk Gregorius

Share on other sites
2 hours ago, Dirk Gregorius said:

You don't need to test edge directions for two AABB. This is only necessary for OBB.

Thanks for the input.

2 hours ago, Dirk Gregorius said:

Then you can clip one against the side planes of the other. This gives you the contact points.

How do I do that? Can you give me a pseudo code for that? It took me several hours to code that SAT part. At the moment, I use this code for finding contact point and normal. Is that right?

static bool TestAxis(XMVECTOR axis, AABB& a, AABB& b, XMVECTOR& mtd, float& mtd2)
{
// intervals along axis
float mina, maxa;
float minb, maxb;
a.supportInterval(axis, mina, maxa);
b.supportInterval(axis, minb, maxb);

// calculate the two possible overlap ranges.
// either we overlap on the left or right sides.
float d0 = (maxb - mina); // 'left' side
float d1 = (maxa - minb); // 'right' side

// the axis length squared
float axis_length_squared = XMVector3Dot(axis, axis).m128_f32[0];

// axis is degenerate. ignore.
if(axis_length_squared < 1.0e-8f)
return true;

// intervals do not overlap. no intersection.
if(d0 < 0.0f || d1 < 0.0f)
return false;

// find out if we overlap on the 'right' or 'left' of the object.
float overlap = (d0 < d1)? d0 : -d1;

// the mtd vector for that axis
XMVECTOR sep = axis * (overlap / axis_length_squared);

// the mtd vector length squared.
float sep_length_squared = XMVector3Dot(sep, sep).m128_f32[0];

// if that vector is smaller than our computed MTD (or the mtd hasn't been computed yet)
// use that vector as our current mtd.
if(sep_length_squared < mtd2 || mtd2 < 0.0f)
{
mtd2 = sep_length_squared;
mtd  = sep;
}
return true;
}

static bool AABBAABBIntersection(RigidBody* rba, AABB& a, RigidBody* rbb, AABB& b, Contact* c)
{
XMVECTOR mtd = XMVectorSet(0, 0, 0, 0), mtddirection;
float mtd2 = -1.0f;

// test faces of A
for(int i=0; i<a.supportFaceCount(); i++)
{
if(!TestAxis(a.supportFaceDirection(i), a, b, mtd, mtd2))
return false;
}

// test faces of B
for(int i=0; i<b.supportFaceCount(); i++)
{
if(!TestAxis(b.supportFaceDirection(i), a, b, mtd, mtd2))
return false;
}

// no intersections were ever tested.
if(mtd2 < 0.0f)
return false;

float Penetration = sqrt(mtd2);

c->a = rba;
c->b = rbb;
c->n = XMVector3Normalize(mtd);
c->p = rba->GetPosition() + mtd*Penetration;

return true;
}

Edited by isu diss

Share on other sites
6 hours ago, isu diss said:

Thanks for the input.

How do I do that? Can you give me a pseudo code for that? It took me several hours to code that SAT part. At the moment, I use this code for finding contact point and normal. Is that right?

static bool TestAxis(XMVECTOR axis, AABB& a, AABB& b, XMVECTOR& mtd, float& mtd2)
{
// intervals along axis
float mina, maxa;
float minb, maxb;
a.supportInterval(axis, mina, maxa);
b.supportInterval(axis, minb, maxb);

// calculate the two possible overlap ranges.
// either we overlap on the left or right sides.
float d0 = (maxb - mina); // 'left' side
float d1 = (maxa - minb); // 'right' side

// the axis length squared
float axis_length_squared = XMVector3Dot(axis, axis).m128_f32[0];

// axis is degenerate. ignore.
if(axis_length_squared < 1.0e-8f)
return true;

// intervals do not overlap. no intersection.
if(d0 < 0.0f || d1 < 0.0f)
return false;

// find out if we overlap on the 'right' or 'left' of the object.
float overlap = (d0 < d1)? d0 : -d1;

// the mtd vector for that axis
XMVECTOR sep = axis * (overlap / axis_length_squared);

// the mtd vector length squared.
float sep_length_squared = XMVector3Dot(sep, sep).m128_f32[0];

// if that vector is smaller than our computed MTD (or the mtd hasn't been computed yet)
// use that vector as our current mtd.
if(sep_length_squared < mtd2 || mtd2 < 0.0f)
{
mtd2 = sep_length_squared;
mtd  = sep;
}
return true;
}

static bool AABBAABBIntersection(RigidBody* rba, AABB& a, RigidBody* rbb, AABB& b, Contact* c)
{
XMVECTOR mtd = XMVectorSet(0, 0, 0, 0), mtddirection;
float mtd2 = -1.0f;

// test faces of A
for(int i=0; i<a.supportFaceCount(); i++)
{
if(!TestAxis(a.supportFaceDirection(i), a, b, mtd, mtd2))
return false;
}

// test faces of B
for(int i=0; i<b.supportFaceCount(); i++)
{
if(!TestAxis(b.supportFaceDirection(i), a, b, mtd, mtd2))
return false;
}

// no intersections were ever tested.
if(mtd2 < 0.0f)
return false;

float Penetration = sqrt(mtd2);

c->a = rba;
c->b = rbb;
c->n = XMVector3Normalize(mtd);
c->p = rba->GetPosition() + mtd*Penetration;

return true;
}

Well, your function computes the minimum translation vector - which is nothing more than a vector you can apply, to resolve a penetration in the most basic way. Dont know what you are making, but if this is enough for you - then its totally fine. For example, for a simple game this may be totally fine. If not, check out dirks great post about contact generation or other resources regarding this topic.

Also, if you want to apply rotation to your bodies, then this system does not work at all - then you definitly need to deal with OBB´s instead.

Share on other sites

@Finalspace I'm making a cricket game. So far, I've been able to resolve collision for ball-ground(ball without any rotation after hitting the ground) & ball-stump(not accurately). I need accurate contact point (c->p) to compute rotation part in the collision resolving code. As you can see the above posted code (c->p) is computed using (c->n), cross product of ra n rb with n goes to zero in the collision resolve code.

void CollisionResponse(Contact *c, float epsilon)
{
XMVECTOR pbdot = c->b->GetVelocityAtPoint(c->p);
XMVECTOR n = c->n;
XMVECTOR ra = (c->p - c->a->GetPosition());
XMVECTOR rb = (c->p - c->b->GetPosition());
XMVECTOR vrel = XMVector3Dot(c->n, (padot - pbdot));
float numerator = (-(1.0f + epsilon)*vrel.m128_f32[0]);
float term1 = (1.0f / c->a->GetMass());
float term2 = (1.0f / c->b->GetMass());
XMVECTOR term3 = XMVector3Dot(c->n, XMVector3Cross(XMVector4Transform(XMVector3Cross(ra, n), c->a->GetIInverse()), ra));
XMVECTOR term4 = XMVector3Dot(c->n, XMVector3Cross(XMVector4Transform(XMVector3Cross(rb, n), c->b->GetIInverse()), rb));
float j = (numerator / (term1+ term2 + term3.m128_f32[0] + term4.m128_f32[0]));
XMVECTOR f = (j*n);

}

Share on other sites
3 hours ago, isu diss said:

@Finalspace I'm making a cricket game. So far, I've been able to resolve collision for ball-ground(ball without any rotation after hitting the ground) & ball-stump(not accurately). I need accurate contact point (c->p) to compute rotation part in the collision resolving code. As you can see the above posted code (c->p) is computed using (c->n), cross product of ra n rb with n goes to zero in the collision resolve code.

void CollisionResponse(Contact *c, float epsilon)
{
XMVECTOR pbdot = c->b->GetVelocityAtPoint(c->p);
XMVECTOR n = c->n;
XMVECTOR ra = (c->p - c->a->GetPosition());
XMVECTOR rb = (c->p - c->b->GetPosition());
XMVECTOR vrel = XMVector3Dot(c->n, (padot - pbdot));
float numerator = (-(1.0f + epsilon)*vrel.m128_f32[0]);
float term1 = (1.0f / c->a->GetMass());
float term2 = (1.0f / c->b->GetMass());
XMVECTOR term3 = XMVector3Dot(c->n, XMVector3Cross(XMVector4Transform(XMVector3Cross(ra, n), c->a->GetIInverse()), ra));
XMVECTOR term4 = XMVector3Dot(c->n, XMVector3Cross(XMVector4Transform(XMVector3Cross(rb, n), c->b->GetIInverse()), rb));
float j = (numerator / (term1+ term2 + term3.m128_f32[0] + term4.m128_f32[0]));
XMVECTOR f = (j*n);

}

Okay then you have no choice and need to produce proper contact points and switch to OBB instead.

For sphere vs OBB you need just one contact point, for OBB vs OBB you can have at max 4 contact points.

I cannot explain how the generation in actual 3D works, but the concept for 2D are similar.

First of, you use SAT to determine the normal with most overlap and record the penetration depth and the normal/axis -> This is what you actually have right now. Also you need to record the body and the face where the axis was found, this is either A or B.

Treat A as the reference and B as the incident body first. But if the most penetrating axis is found on B, then B is your reference body and A is your incident body. Its just swapped.

The recorded body and face A is your reference face. The most anti-parallel face on B is your incident face.

Next is to clip the incident face against the side planes of the reference face (Sutherland Hodgman clipping) -> You dont clip against the reference face.

You save all the clipped points below the reference face. Note that this clipped points are on the incident face - not on the reference face!

For better coherence you project the saved clip points to the reference face.

Now you have all the contact points on A. To get the second contact point its as simple as: PointOnB = PointOnA + -Normal * PenetrationDistance;

There you have it.

Again, everything i am talking about is good explained in "DirkGregorius_Contacts.pdf" - even with code!

One last tip i can give you: Visualize everything! Draw the face planes, the clipping planes, the contact points and the penetration depth. Colorize the reference body red, the incident body blue  Trust me, this helps a lot!

Edited by Finalspace

Share on other sites

I learned SAT and applied for AABB. It gives me part of the answer. To have the full answer, I need S-H clippings. But S-H clipping is so hard. @Finalspacetold me how to combine SAT result with S-H clippings for AABB. I don't have much experience of implementing like that. Can anybody share sourcecode? even one bit? please, I wouldn't ask if I knew

Share on other sites

Here you go:

//--------------------------------------------------------------------------------------------------
RnArray< RnVector3 > rnClipPolygon( const RnArray< RnVector3 >& Polygon, const RnPlane& Plane, int Edge )
{
RN_ASSERT( Polygon.Size() >= 3 );
RnArray< RnVector3 > Out;

RnVector3 Vertex1 = Polygon.Back();
float Distance1 = rnDistance( Plane, Vertex1 );

for ( int Index = 0; Index < Polygon.Size(); ++Index )
{
RnVector3 Vertex2 = Polygon[ Index ];
float Distance2 = rnDistance( Plane, Vertex2 );

if ( Distance1 <= 0.0f && Distance2 <= 0.0f )
{
// Both vertices are behind the plane - keep vertex2
Out.PushBack( Vertex2 );
}
else if ( Distance1 <= 0.0f && Distance2 > 0.0f )
{
// Vertex1 is behind of the plane, vertex2 is in front -> intersection point
float Fraction = Distance1 / ( Distance1 - Distance2 );
RnVector3 IntersectionPoint = Vertex1.Position + Fraction * ( Vertex2.Position - Vertex1.Position );

// Keep intersection point
Out.PushBack( IntersectionPoint );
}
else if ( Distance2 <= 0.0f && Distance1 > 0 )
{
// Vertex2 is behind of the plane, vertex1 is in front -> intersection point
float Fraction = Distance1 / ( Distance1 - Distance2 );
RnVector3 IntersectionPoint = Vertex1.Position + Fraction * ( Vertex2.Position - Vertex1.Position );

// Keep intersection point
Out.PushBack( IntersectionPoint );

// And also keep vertex2
Out.PushBack( Vertex2 );
}

// Keep vertex2 as starting vertex for next edge
Vertex1 = Vertex2;
Distance1 = Distance2;
}

return Out;
}

HTH,

-Dirk

Create an account

Register a new account

• What is your GameDev Story?

In 2019 we are celebrating 20 years of GameDev.net! Share your GameDev Story with us.

• 11
• 13
• 9
• 25
• 18
• Forum Statistics

• Total Topics
634437
• Total Posts
3017459
×


text:
# Post-processed quantities

## Introduction

This section covers quantities, which can be evaluated by means of the Flux postprocessor. These quantities are:

• local (or spatial ) quantities, which can be evaluated in each point of the computation domain
• global quantities, which are results of integration of local quantities over the entire or selected regions of the computation domain

format: plain


text:
# Completing the square

Completing the square is a way to solve a quadratic equation if the equation will not factorise.

It is often convenient to write an algebraic expression as a square plus another term. The other term is found by dividing the coefficient of $$x$$ by $$2$$, and squaring it.

Any quadratic equation can be rearranged so that it can be solved in this way.

When completing the square, we end up with the form:

$y=a(x\pm b)^2+c$


text:
start:development:ides_and_editors

# Differences

This shows you the differences between two versions of the page.

 start:development:ides_and_editors [2021/10/12 11:27]jrutte02 start:development:ides_and_editors [2021/10/12 11:36] (current)jrutte02 2021/10/12 11:36 jrutte02 2021/10/12 11:27 jrutte02 2020/04/16 19:29 jrutte02 created 2021/10/12 11:36 jrutte02 2021/10/12 11:27 jrutte02 2020/04/16 19:29 jrutte02 created Line 1: Line 1: - ====== Editors and IDEs ====== + - + - ===== Editors ===== + - + ===== Editors, IDEs and Notebooks ===== - + - On our clusters we support a number of editors. If you think your favorite one is missing, [[https://hpc.uni-mainz.de/high-performance-computing/service-angebot/softwareinstallation/|please let us know]]. + - + - The following list is not exhaustive nor complete. + - + * {{fa>file-code-o?fw&color=#000}} [[:start:development:ides_and_editors:editors|Editors]] - + * {{fa>magic?fw&color=#000}} [[:start:development:ides_and_editors:ides|IDEs]] - + * {{fa>line-chart?fw&color=#000}} [[:start:development:ides_and_editors:notebooks|Notebooks]] - ==== atom ==== + - + - [[https://atom.io/|''atom'']] is a widely used IDE, favored by many developers, which supports all the widespread programming languages and provides a plug-in scheme for most purposes. + - + - ''atom'' can be installed on various platforms. Please refer to the [[https://flight-manual.atom.io/getting-started/sections/installing-atom/|installation guide]] (see the top tap for platforms other than Linux). + - + - === Running Locally and Editing Remote Files === + - + - If you would like to develop on our HPC-Systems, you can edit the source code files stored there. We recommend using the ''ftp-remote-edit'' package. Select: ''Edit'' -> ''Preferences'' -> ''Install'' to select this package. + - + - After installing ''ftp-remote-edit'' you can select it in the ''Packages'' drop-down menue and select ''Toggle''. You will be asked for a password, which is actually a passphrase you may choose at will. + - + - You will need to enter this password //only once per session//, much like for ordinary [[:start:mogon_cluster:access|ssh-connections]]. + - + - + - + - == Working within the University Network == + - + - After the installation of ''ftp-remote-edit'' has been successfully completed, the following settings have to be made in order to work within the University Network. + - Select: ''Packages'' -> ''Ftp-remote-edit'' -> ''Toggle'' and enter your password. On the left side of the editor the remote pane has now opened. Right click on the pane and select ''Edit Servers''. Enter the settings for the desired cluster according the following template: + - + - ^ MOGON I ^ MOGON II ^ + - | {{ :development:atom_ftp_remote_edit_mi.png?400 | example setting for MOGON I }}  | {{ :development:atom_ftp_remote_edit_mii.png?400 | example setting for MOGON II}}  | + - + - + - == Working Remotely outside the University Network == + - + - Thereafter, selecting ''Toggle'' will give you the remote pane with a ''Edit Servers'' button. Select it and enter settings according to this template: + - + - {{:development:atom.png?400 | example setting for the MOGON I fileserver}} + - + - Please note: + - * Select the hostname for direct connect to the file server. This can be either ''mogonfs.zdv.uni-mainz.de'' for MOGON I or ''mogon2ftp.zdv.uni-mainz.de''. If you develop in your home directory, this choice does not matter. See the [[:start:fs_dm:filesystems:transfer|filesystem page]] for further details. + - * Choose the plain FTP protocol, but require an encrypted handshake with ''Require explicit FTP over TLS''. + - * As the logon type should be ''Username / Password'', supply them in the according fields. + - * Finally supply a full qualified path, e.g. ''/gpfs/fs1/home//...'' or ''/lustre/miifs01/project/''. See the [[:start:fs_dm:filesystems:transfer|filesystem page]] for further options. + - + - + - + - ==== gedit ==== + - + - [[https://wiki.gnome.org/Apps/Gedit|gedit]] is a nice and simple to use editor with syntax highlighting for a number of languages. Due to its simplicity it is the favorite tool for some of our courses. + - + - === Turning off Gtk Warnings === + - + - gedit issues a number of (annoying and meaningless) Gtk warnings. They can be turned off by placing this snippet in your ''.bashrc''. + - + - + - suppress-gnome-warnings() { + - # $1 is the name which should appear on history but is otherwise unused. + - historyName=$1 + - shift + - + - if [ -n "$*" ]; then + - # write the real command to history without the prefix + - history -s "$historyName ${@:2}" + - + - # catch the command output + - errorMsg=$( $* 2>&1 ) + - + - # check if the command output contains not a (one of two) GTK-Warnings + - if !$(echo $errorMsg | grep -q 'Gtk-WARNING\|connect to accessibility bus'); then + - echo$errorMsg + - fi + - fi + - } + - + - gedit() { + - suppress-gnome-warnings $FUNCNAME$(which $FUNCNAME)$@ + - } + - + - + - + - ==== emacs ==== + - + - ''emacs'' is provided on every login node. + - + - + - + - ==== nano ==== + - + - ''nano'' is provided on every login node. + - + - + - ==== vi / vim ==== + - + - ''vi'' (as ''vim'') is provided on every login node. + - + - ===== IDEs ===== + - + - ==== Geany ==== + - + - [[https://www.geany.org/|''geany'']] is provided a system package. You can simply use it on any login node. + - + - ====== Notebooks ====== + - + - * [[start:development:ides_and_editors:notebooks|Notebooks]] +
• start/development/ides_and_editors.txt


text:
## Coarsening dynamics of three-dimensional levitated foams : From wet to dry

2013
Journal article
##### Published in
The European Physical Journal E ; 36 (2013), 10. - 116. - ISSN 1292-8941. - eISSN 1292-895X
##### Abstract
We study diamagnetically levitated foams with widely different liquid fractions. Due to the levitation, drainage is effectively suppressed and the dynamics is driven by the coarsening of the foam bubbles. For dry foams, the bubble size is found to increases as the square root of foam age, as expected from a generalized von Neumann law. At higher liquid content the behavior changes to that of Ostwald ripening where the bubbles grow with the 1/3 power of the age. Using Diffusing Wave Spectroscopy we study the local dynamics in the different regimes and find diffusive behavior for dry foams and kinetic behavior for wet foams.
530 Physics
##### Cite This
ISO 690ISERT, Nathan, Georg MARET, Christof M. AEGERTER, 2013. Coarsening dynamics of three-dimensional levitated foams : From wet to dry. In: The European Physical Journal E. 36(10), 116. ISSN 1292-8941. eISSN 1292-895X. Available under: doi: 10.1140/epje/i2013-13116-x
BibTex
@article{Isert2013-10Coars-26658,
year={2013},
doi={10.1140/epje/i2013-13116-x},
title={Coarsening dynamics of three-dimensional levitated foams : From wet to dry},
number={10},
volume={36},
issn={1292-8941},
journal={The European Physical Journal E},
author={Isert, Nathan and Maret, Georg and Aegerter, Christof M.},
note={Article Number: 116}
}

RDF
<rdf:RDF
xmlns:dcterms="http://purl.org/dc/terms/"
xmlns:dc="http://purl.org/dc/elements/1.1/"
xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
xmlns:bibo="http://purl.org/ontology/bibo/"
xmlns:dspace="http://digital-repositories.org/ontologies/dspace/0.1.0#"
xmlns:foaf="http://xmlns.com/foaf/0.1/"
xmlns:void="http://rdfs.org/ns/void#"
xmlns:xsd="http://www.w3.org/2001/XMLSchema#" >
<dcterms:available rdf:datatype="http://www.w3.org/2001/XMLSchema#dateTime">2014-03-04T15:29:56Z</dcterms:available>
<dcterms:hasPart rdf:resource="https://kops.uni-konstanz.de/bitstream/123456789/26658/2/Isert_266585.pdf"/>
<dcterms:abstract xml:lang="eng">We study diamagnetically levitated foams with widely different liquid fractions. Due to the levitation, drainage is effectively suppressed and the dynamics is driven by the coarsening of the foam bubbles. For dry foams, the bubble size is found to increases as the square root of foam age, as expected from a generalized von Neumann law. At higher liquid content the behavior changes to that of Ostwald ripening where the bubbles grow with the 1/3 power of the age. Using Diffusing Wave Spectroscopy we study the local dynamics in the different regimes and find diffusive behavior for dry foams and kinetic behavior for wet foams.</dcterms:abstract>
<dc:creator>Maret, Georg</dc:creator>
<dc:contributor>Isert, Nathan</dc:contributor>
<dcterms:issued>2013-10</dcterms:issued>
<foaf:homepage rdf:resource="http://localhost:8080/"/>
<dspace:isPartOfCollection rdf:resource="https://kops.uni-konstanz.de/server/rdf/resource/123456789/41"/>
<dc:creator>Isert, Nathan</dc:creator>
<dcterms:rights rdf:resource="https://rightsstatements.org/page/InC/1.0/"/>
<dspace:hasBitstream rdf:resource="https://kops.uni-konstanz.de/bitstream/123456789/26658/2/Isert_266585.pdf"/>
<dcterms:isPartOf rdf:resource="https://kops.uni-konstanz.de/server/rdf/resource/123456789/41"/>
<dc:contributor>Maret, Georg</dc:contributor>
<dcterms:bibliographicCitation>European Physical Journal E ; 36 (2013), 10. - 116</dcterms:bibliographicCitation>
<dcterms:title>Coarsening dynamics of three-dimensional levitated foams : From wet to dry</dcterms:title>
<dc:rights>terms-of-use</dc:rights>
<dc:date rdf:datatype="http://www.w3.org/2001/XMLSchema#dateTime">2014-03-04T15:29:56Z</dc:date>
<void:sparqlEndpoint rdf:resource="http://localhost/fuseki/dspace/sparql"/>
<dc:contributor>Aegerter, Christof M.</dc:contributor>
<bibo:uri rdf:resource="http://kops.uni-konstanz.de/handle/123456789/26658"/>
<dc:language>eng</dc:language>
<dc:creator>Aegerter, Christof M.</dc:creator>
</rdf:Description>
</rdf:RDF>

Yes


text:
# Transient Reduction Pipeline

2010 June 20 - T. Joseph W. Lazio

There are a class of observations for which only relatively simple data reduction steps are needed. One such example is that of transient observations, which are typically conducted in continuum mode and for which one is merely trying to determine a flux density of an essentially unresolved or only partially resolved source. This guide describes the steps in such a pipeline.

In order to process these data in a semi-automatic fashion, certain assumptions are made

• Data stored as an Archive Science Data Model (ASDM) file on disk.
• Sources, listed in order of appearance, are structured as Option A or Option B
• Option A: phase calibrator, target source, amplitude calibrator
• Option B: amplitude calibrator, phase calibrator, target source

This guide should be set up so that a pipeline script can be extracted from it.

Warning: This guide was written at the time when the EVLA was still in its commissioning phase. As the instrument matures, specific steps taken here may need to be adjusted. Caveat emptor. Possible issues include

• calwt=F because at the time of writing the EVLA is not reporting accurate weights (though I believe that this has a stronger effect on extended sources rather than compact sources);
• the pipeline assumes the same number of channels per spectral window;
• the pipeline assumes that polarization information is not required;
• the pipeline assumes that the first scan is a dummy scan, as such is required at the time of writing.

Further, to the best of the author's knowledge, the infrastructure does not yet exist within CASA to determine if there are antenna position corrections that should be applied (via gencal).

I've restructured the document to describe the processing steps (conceptually) first, then list the actual script itself.

## User Input

I am attempting to migrate the script to become a "pipeline," in the sense of not requiring any user input. However, there are some pieces of information that it is useful or required to know in order to process the data.

• What is the name of the initial Archive Science Data Model (ASDM) file? It would be possible to assume that the script is being run in a directory that contains only a single file, which is in ASDM format, but that also seems a bit limiting.
• What antenna should be used as a reference antenna? At the time of writing, it is not clear that a robust algorithm exists within CASA for choosing a reference antenna if the user has not specified one.
• How much flagging should be done? The script does some simple clipping, designed to remove any egregious RFI or horribly performing antenna or baseline. This clipping is done in terms of the rms visibility amplitude in the science target data, flagging data above some large threshold (e.g., 50$\sigma$). This threshold is under user control, but the current flagging in this script is certainly not equal to a human lovingly massaging the visibility data.

These must be done before invoking the script/pipeline. For example

CASA <1> ASDM_name='AL007_sb123456789_1.56000.1234567890'
CASA <2> refant='ea21'
CASA <3> clip_sigma=50.0
CASA <4> execfile('transient_pipeline.py')

A sensible, large value is adopted for clip_sigma if it is omitted.

## Processing Steps

### Initial Stuff

• Read the data from the ASDM file converting it to a Measurement Set with importevla. Apply basic flagging operations (zeros, shadowing) here.
• Flag the first (dummy) scan. (This step is required at the time of development, but it may be relaxed in the future.)
• Flag (quack) first 10 seconds of each scan.
• Having constructed the initial measurement set, pause to extract various useful items from it, such as the frequency of observation, number of spectral channels, etc. Use these to then calculate various quantities such as the primary beam, synthesized beam, tolerable amount of bandwidth smearing, and tolerable amount of time-average smearing.
• Reduce the data in size for faster processing downstream by averaging in time and frequency. Also, reject the edge channels.

### Calibration

• Set the flux density of the amplitude calibrator, if one exists. The pipeline uses the calibration codes stored in the measurement set to determine which, if any, sources are suitable for calibrating the flux density scale. If it finds no source that is suitable, it sets all sources that are suitable for calibrating the gain phases and the spectral bandpass to have a flux density of 1 Jy. There is a potential issue, at the time of writing, if the frequency setting is such that one should use a model image at a different band than one is observing, e.g., observing near the top of the C band where an X band model might be more appropriate. Presumably once all of the new receivers have been installed, new calibrator models will follow at some point.
• Make a quick-n-dirty bandpass, if a spectral bandpass calibrator exists.
• Amplitude and phase calibration. This is done in two steps, first the amplitude calibrator, then the phase calibrator, but the two steps could be combined.
• Apply the calibration, from the phase calibrator to the science target.
• Form target source measurement set.
• Clip data above some a threshold. The user can select this threshold by specifying clip_sigma prior to invoking the pipeline. If no value is given, a value of 50$\sigma$ is assumed.

### Imaging

• The source is imaged and a light CLEANing is applied.

## Pipeline

# In CASA

#!/bin/env python
#
# A first step toward an EVLA transient reduction pipeline.
#
# Warning: This script was written at the time that the
# EVLA was undergoing commissioning.  Details may change.
# Caveat emptor!
#
# In particular, current version does not calibrate
# weights within applycal as EVLA does not currently
# report accurate weights.
#
#
########################################################################
#
import os.path, math
#
#Version='v. 0.0 TJWL  2010-05-11'
#Version='v. 0.1 TJWL  2010-06-16'
#Version='v. 0.2 TJWL  2010-06-19'
Version='v. 0.3 TJWL  2010-06-20'
# adopt L. Chiomuk's suggestion of incorporating spectral window
# information
# start moving away from user input
#
CalModels='/home/casa/data/nrao/VLA/CalModels'
c=2.99792458E8
D=24.5
arcsecond_deg=3600
#
print " "
print "EVLA Transient Reduction Pipeline"
print " "
print Version
print " "
print "Warning: This script was written at the time that the"
print "EVLA was undergoing commissioning.  Details may change."
print "Caveat emptor!"
print " "
print "Assumption: Data stored as an Archive Science Data Model (ASDM)"
print "            file on disk."
print " "
print "Assumption: Sources, listed in order of appearance, are structured"
print "            as Option A or Option B"
print "            Option A: phase calibrator"
print "                      target source"
print "                      amplitude calibrator"
print " "
print "            Option B: amplitude calibrator"
print "                      phase calibrator"
print "                      target source"
print " "
#
print "Assumption: Calibrator model files stored in "
print CalModels
if not os.path.isdir(CalModels):
raise IOError('Calibrator model files missing?  Stopping.')
#

########################################################################
# User specified parameters
try: clip_sigma
except NameError:
clip_sigma=50.0

try: ASDM_name
except NameError:
raise NameError('ASDM_name variable not specified')

if not os.path.exists(ASDM_name):
raise IOError('ASDM file does not exist.')
else:
msname=ASDM_name+'.ms'

try: refant
except NameError:
raise NameError('Reference antenna not specified in variable refant.')

########################################################################
# definitions
#
acal_std_name = {}
acal_std_name['3C286'] = frozenset(['3C286','3C 286',
'1328+30','1328+307','B1328+307','B1328+30',
'1331+305','J1331+305','J1331+3030'])
acal_std_name['3C138'] = frozenset(['3C138','3C 138',
'0518+16','0518+165','0518+1635','B0518+16','B0518+165','B0518+1635',
'0521+166','0521+1638','J0521+166','J0521+1638'])
acal_std_name['3C147'] = frozenset(['3C147','3C 147',
'0538+49','0538+498','0538+4949','B0538+49','B0538+498','B0538+4949',
'0542+498','0542+4951','J0542+498','J0542+4951'])
acal_std_name['3C48'] = frozenset(['3C48','3C 48',
'0134+32','0134+329','0134+3254','B0134+32','B0134+329','B0134+3254',
'0137+331','0137+3309','J0137+331','J0137+3309'])

########################################################################
#
# 2. Determine the structure of the observations.
#
print " "
obs_structure=raw_input('Structure of observations [A|B]: ')
if obs_structure == 'A':
sources = {'acal':'2', 'pcal':'0', 'target':'1'}
elif obs_structure == 'B':
sources = {'acal':'0', 'pcal':'1', 'target':'2'}

#
########################################################################
#
# Initial stuff.
#
# Read the data from the ASDM file converting it to a Measurement Set
# with importevla.
# Apply basic flagging operations (zeros, shadowing) here.
#
importevla(asdm=ASDM_name, vis=msname,
#
# Flag the first (dummy) scan.
# Required at the time of development, this step may be relaxed
# in the future.
#
print "Flagging first (dummy) scan ..."
flagdata(vis=msname,
mode='manualflag',selectdata=T, scan='1')

#
# Flag (quack) first 10 seconds of each scan.
#
print "Quack data ..."
flagdata(vis=msname,
mode='quack',quackinterval=10,quackmode='beg',
selectdata=F)

#
# Extract various useful quantities ....
#
hdkey='field', hdindex=sources['acal'])
acal_name=hdvalue[0]

spw_table=msname + '/SPECTRAL_WINDOW'
tb.open(spw_table)
freq_list=tb.getcol("REF_FREQUENCY")
channel_width=tb.getcol("CHAN_WIDTH")
num_chan=tb.getcol("NUM_CHAN")
tb.close(spw_table)

freq=freq_list[0]
if (1E9 < freq) and (freq < 2E9):
band = 'L'
elif (2E9 < freq) and (freq < 4E9):
band = 'S'
elif (4E9 < freq) and (freq < 8E9):
band='C'
elif (8E9 < freq) and (freq < 12E9):
band='X'
elif (12E9 < freq) and (freq < 40E9):
band = 'K'
elif freq > 40E9:
band = 'Q'
print "Observations are determined to be in the ", band, " band.\n"

wavelength=c/freq
primary_beam=(1.02*wavelength/D)*degrad   # HPBW, from Napier (1999)
FoV=(primary_beam/math.sqrt(2))*arcsecond_deg

tb.open(msname)
uvw=tb.getcol('UVW')
interval=tb.getcol('INTERVAL')
tb.close(msname)

b=[]
for u, v, w in zip(uvw[0], uvw[1], uvw[2]):
b.append(math.sqrt(u*u+v*v)/wavelength)
bmax=max(b)
# HPBW, synthesized beam,
# from Bridle & Schwab (1999)

#
# How much time-average smearing can be tolerated?
# Assume no more than time_smearing_loss peak
# intensity loss over half of primary beam.
# Follow Section 2 of Bridle & Schwab (1999).
# Note undocumented option that the amount of time-average
# smearing loss can be specified by user.
#
try: time_smearing_loss
except NameError:
time_smearing_loss=0.01

tau=math.sqrt(time_smearing_loss*1E9)*(synthesized_beam/FoV)

arbitrary_maximum=30
if tau > arbitrary_maximum:
tau=arbitrary_maximum

dt=min(interval)
#
# tau is allowed value, dt is actual (minimum)
# (could be an issue if baseline-dependent
#  correlator accumulation used)
# make sure that tau is an integer multiple of dt
#
tau=dt*math.floor(tau/dt)
if tau < dt:
tau=dt

print "Data will be averaged in time."
print "Original time sampling [s]: ", dt
print "Averaging time [s]:         ", tau
print " "

#
# How much bandwidth smearing can be tolerated?
# Assume no more than band_smearing_loss peak
# intensity loss over half of primary beam.
# Follow Section 1 of Bridle & Schwab (1999).
# Assume square bandpass, no taper, expand resulting
# sine integral to lowest order
# Note undocumented option that the amount of time-average
# smearing loss can be specified by user.
#
try: band_smearing_loss
except NameError:
band_smearing_loss=0.01

eta_band=3.79
delta_nu=freq*(2/eta_band)*(synthesized_beam/FoV)*math.sqrt(18*band_smearing_loss)
#
# delta_nu is allowed value,
# figure out even divisibles of the actual
# value, stored in num_chan that is smaller than
# delta_nu
# potential bug if uneven channel widths
# used in different spectral windows

dnu=channel_width[0][0]
nchan_log2=math.log(num_chan[0],2)
for i in range(int(nchan_log2-1), 0, -1):
if (dnu*math.pow(2,i)) < delta_nu:
nchav=math.pow(2,i)
break

if nchav < 1:
nchav=1

print "Data will be averaged in frequency."
print "Original channel width [kHz]: ", min(channel_width[0])/1E3
print "Averaged channel width [kHz]: ", (nchav*min(channel_width[0]))/1E3
print "Number of channels: ", nchav

#
# Compress data for faster processing.
#    Throw away edge channels.
#

bchan=num_chan[0]*0.1
echan=num_chan[0]*0.9
chans='*:%d~%d'%(bchan, echan)

outmsname=ASDM_name+'_split.ms'
tave='%.fs'%tau
split(vis=msname,outputvis=outmsname,
datacolumn='data',
spw=chans,
timebin=tave,width=int(nchav))

########################################################################
#
# Calibration
#
# Set flux density of amplitude calibrator.
#    Need some heuristics here.
#    Also could be an issue if the frequency
#    setting is such that one should use a
#    model image at a different band than one
#    is observing, e.g., observing near the
#    top of the C band where an X band model
#    might be more appropriate.
#
#print acal_name
print " Beginning calibration ..."
clearstat()

for tstname in acal_std_name.keys():
if acal_name in acal_std_name[tstname]:
modimage=tstname
modimage=CalModels + '/' + modimage + '_' + band + '.im'
setjy(vis=outmsname,field=sources['acal'],modimage=modimage)

#
# Quick-n-dirty bandpass
#
print " Finding bandpass ..."

Bcaltable=ASDM_name + '.B1'
Gcaltable=ASDM_name + '.G0'
gaincal(vis=outmsname,caltable=Gcaltable,
field=sources['acal'],
solint='int',refant=refant,
solnorm=T,calmode='p',minsnr=5)
bandpass(vis=outmsname,caltable=Bcaltable,
field=sources['acal'],
solint='inf',combine='scan',refant=refant,
solnorm=T,
gaintable=[Gcaltable])

#
# Amplitude and phase calibration
#
print " Amplitude and phase calibration ..."

Gcaltable=ASDM_name + '.G1'
gaincal(vis=outmsname,caltable=Gcaltable,
field=sources['acal'],
solint='inf',combine='scan',refant=refant,
minsnr=5,
gaintable=[Bcaltable])

gaincal(vis=outmsname,caltable=Gcaltable,
field=sources['pcal'],
solint='inf',refant=refant,
minsnr=5,
append=T,
gaintable=[Bcaltable])

#
# Apply calibration.
#
fluxtable=ASDM_name + '.flux1'
fluxscale(vis=outmsname,caltable=Gcaltable,fluxtable=fluxtable,
reference=sources['acal'],
transfer=sources['pcal'])

applycal(vis=outmsname,
field=sources['pcal'],
gaintable=[Bcaltable, fluxtable],
gainfield=['',sources['pcal']],
calwt=F)

applycal(vis=outmsname,
field=sources['target'],
gaintable=[Bcaltable, fluxtable],
gainfield=['',sources['pcal']],
calwt=F)

#
# Form target source measurement set.
#
print " Forming target source measurement set ..."

hdkey='field', hdindex=sources['target'])
target_name=hdvalue[0]

targetms=target_name + '.ms'
split(vis=outmsname,outputvis=targetms,
datacolumn='corrected',
field=sources['target'])

#
# Clip any egregious data.
#
Vstat=visstat(vis=outmsname,selectdata=F)
Vclip=clip_sigma*Vstat['DATA']['rms']

flagdata(vis=outmsname,
mode='manualflag',
clipexpr='ABS RR',clipminmax=[0.0,Vclip],
selectdata=F)
flagdata(vis=outmsname,
mode='manualflag',
clipexpr='ABS LL',clipminmax=[0.0,Vclip],
selectdata=F)

########################################################################
#
# Imaging and CLEANing.
#
print " Image the data ..."
clearstat()

imagename=target_name + '.im'
cells='%.3farcsec'%(synthesized_beam/4)

clean(vis=targetms,imagename=imagename,
mode='mfs',
niter=300,
psfmode='clark',
imagermode='csclean',
imsize=[1024,1024],cell=cells,
stokes='IV',
weighting='briggs',robust=-1)


text:
## Synopsis

The ZIPload model is designed to provide a means to access the Residential_enduse class at a fundamental level and from a "power engineer" perspective. It uses a classic ZIP load model (constant impedance, current, and power) where "base power" is specified, then the ZIP fractions and power factors are assigned. The ZIP model equations can be defined as:

$P_i = \frac{|V_a^2|}{|V_n^2|}*|S_n|*Z_{\%}*cos(Z_{\theta}) + \frac{|V_a|}{|V_n|}*|S_n|*I_{\%}*cos(I_{\theta}) + |S_n|*P_{\%}*cos(P_{\theta})$

$Q_i = \frac{|V_a^2|}{|V_n^2|}*|S_n|*Z_{\%}*sin(Z_{\theta}) + \frac{|V_a|}{|V_n|}*|S_n|*I_{\%}*sin(I_{\theta}) + |S_n|*P_{\%}*sin(P_{\theta})$

where:

• $P_i$: Real power consumption of the ith load
• $Q_i$: Reactive power consumption of the ith load
• $V_a$: Actual terminal voltage
• $V_n$: Nominal terminal voltage
• $S_n$: Apparent Power consumption at nominal voltage
• $Z_\%$: Percent of load that is constant impedance
• $I_\%$: Percent of load that is constant current
• $P_\%$: Percent of load that is constant power
• $Z_\theta$: Phase angle of constant impedance fraction
• $I_\theta$: Phase angle of constant current fraction
• $P_\theta$: Phase angle of constant power fraction

In a time-variant load representation, the coefficients of the ZIP model, $V_n, S_n, Z_\%, I_\%, P_\%, Z_\theta, I_\theta,$ and $P_\theta$ remain constant, but the power consumption, $P_i$ and $Q_i$, of the ith load varies with the actual terminal voltage, $V_a$. The ZIP model is similar to the polynomial representation used in many commercial software packages. In the polynomial representation of the ZIP load, the constant coefficient is roughly equivalent to the P fraction, the linear coefficient is roughly equivalent to the I fraction, and the quadratic coefficient is roughly equivalent to the Z fraction. The ZIP model only varies the power consumption as a function of actual terminal voltage, $V_a$. Note also that the sum of $Z_\%, I_\%,$ and $P_\%$ must equal one.

Internally, the model performs all of the proper phase rotations and scaling for various voltage levels (e.g., the model would look the same for a house attached on Phase A versus Phase C). Currently, it is assumed that nominal voltage is either 120 or 240 V (this will be modified in future upgrades to residential_enduse).

## Parameters

Property Name Type Unit Description
base_power double kW Base real power of the total load at nominal voltage.
power_pf double pu Power factor for constant power portion of load.
current_pf double pu Power factor for constant current portion of load.
impedance_pf double pu Power factor for constant impedance portion of load.
actual_power complex kVA [read only] Variable to monitor total power of load as a function of voltage.
heatgain_only boolean - Toggles the zipload to generate heat only (no kW), is deactivated (false) by default
is_240 boolean - Toggles between a 120 V (unbalanced) vs. a 240 V (balanced) connection - true indicates it is a 240 V load.
These variables are inherited from the enduse load structure.
power_fraction double pu The fraction of the load that is constant power.
current_fraction double pu The fraction of the load that is constant current.
impedance_fraction double pu The fraction of the load that is constant impedance.
heatgain_fraction double pu The fraction of the total load (kW) that produces waste heat.
These variables are used for creating cyclic load behavior.
duty_cycle double pu Fraction of time the device is in the on state.
period double hours Time interval over which duty cycle is applied.
phase double pu Indicates percentage of phase that the object is currently in; running period is assumed to be from 0.0 until percent of duty cycle. To create a distribution of devices, this variable should be randomized between 0 and 1.
These variables are used for various aggregate and demand response modes.
demand_response_mode boolean - Activates equilibrium dynamic representation of demand response.
number_of_devices int16 - [Used with demand_response_mode only.] Number of devices to model - base power is the total load of all devices.
thermostatic_control_range int16 K [Used with demand_response_mode only.] Range of the thermostat's control operation.
number_of_devices_off double - [Used with demand_response_mode only.] Total number of devices that are off.
number_of_devices_on double - [Used with demand_response_mode only.] Total number of devices that are on.
rate_of_cooling double K/h [Used with demand_response_mode only.] rate at which devices cool down.
rate_of_heating double K/h [Used with demand_response_mode only.] rate at which devices heat up.
temperature int16 K [Used with demand_response_mode only.] temperature of the device's controlled media (eg air temp or water temp).
phi double pu [Used with demand_response_mode only.] duty cycle of the device(s).
demand_rate double 1/h [Used with demand_response_mode only.] consumer demand rate that prematurely turns on a device or population.
nominal_power double kW [Used with demand_response_mode only.] the rated amount of power demanded by devices that are on.
recovery_duty_cycle double pu [Used with cycling mode and passive_controller duty cycle mode.] Fraction of time in the on state, while in recovery interval.
multiplier double pu [Used with cycling mode and passive_controller duty cycle mode.] This variable is used to modify the base power as a function of multiplier times base_power.

TODO: : Document cycling, demand response and aggregate modes.

## Example models

This model is representative of a simple ZIPload with scheduled load behavior.

    object ZIPload {
parent house1;
heatgain_fraction 0.90;
power_pf 1.0;
current_pf 1.0;
impedance_pf 1.0;
impedance_fraction 0.20;
current_fraction 0.40;
power_fraction 0.40;
};

This model is representative of a ZIPload used in cycling mode to roughly represent a schedule pool pump (note, pool_pump_season is a schedule of 0 or 1).

    object ZIPload {
name house1_poolpump;
parent house1;
base_power pool_pump_season*1.44;
duty_cycle 0.22;
phase 0.26;
period 4.96;
heatgain_fraction 0.0;
power_pf 1.0;
current_pf 1.0;
impedance_pf 1.0;
impedance_fraction 0.20;
current_fraction 0.40;
power_fraction 0.40;
is_240 TRUE;
};

This model is representative of a ZIPload with a passive controller used to implement the elasticity model out of the market module.

    object ZIPload {
parent house1;
heatgain_fraction 0.90;
power_pf 1.0;
current_pf 1.0;
impedance_pf 1.0;
impedance_fraction 0.20;
current_fraction 0.40;
power_fraction 0.40;
object passive_controller {
period 900;
control_mode ELASTICITY_MODEL;
two_tier_cpp true;
observation_object Market_1;
observation_property past_market.clearing_price;
state_property multiplier;
linearize_elasticity true;
price_offset 0.01;
critical_day CPP_days_R1.value;
first_tier_hours 12;
second_tier_hours 12;
third_tier_hours 6;
first_tier_price 0.060483;
second_tier_price 0.120965;
third_tier_price 0.604826;
old_first_tier_price 0.124300;
old_second_tier_price 0.124300;
old_third_tier_price 0.124300;
daily_elasticity daily_elasticity_wtech*1.1731;
sub_elasticity_first_second -0.1783;
sub_elasticity_first_third -0.2604;
};
};

This model is representative of a "pool pump" or cycling model that is using a DR control from the market module.

    object ZIPload {
name house1_poolpump;
parent house1;
base_power pool_pump_season*1.44;
duty_cycle 0.22;
phase 0.26;
period 4.96;
heatgain_fraction 0.0;
power_pf 1.0;
current_pf 1.0;
impedance_pf 1.0;
impedance_fraction 0.20;
current_fraction 0.40;
power_fraction 0.40;
is_240 TRUE;
recovery_duty_cycle 0.27;
object passive_controller {
period 900;
control_mode DUTYCYCLE;
pool_pump_model true;
observation_object Market_1;
observation_property past_market.clearing_price;
state_property override;
base_duty_cycle 0.22;
setpoint duty_cycle;
first_tier_hours 12;
second_tier_hours 12;
third_tier_hours 6;
first_tier_price 0.060483;
second_tier_price 0.120965;
third_tier_price 0.604826;
};
};

TODO: : Examples for cycling, demand response and aggregate modes.


text:
# Greatest Greatest Common Divisor

Time Limit: 4000/2000 MS (Java/Others)

Memory Limit: 65536/65536 K (Java/Others)

## Description

Pick two numbers ${a}_{i},{a}_{j}(i\neq j)$ from a sequence to maximize the value of their greatest common divisor.

## Input

Multiple test cases. In the first line there is an integer $T$, indicating the number of test cases. For each test cases, the first line contains an integer $n$, the size of the sequence. Next line contains $n$ numbers, from ${a}_{1}$ to ${a}_{n}$. $1\leq T\leq 100,2\leq n\leq {10}^{5},1\leq {a}_{i}\leq {10}^{5}$. The case for $n\geq {10}^{4}$ is no more than $10$.

## Output

For each test case, output one line. The output format is Case #$x$: $ans$, $x$ is the case number, starting from $1$, $ans$ is the maximum value of greatest common divisor.

## Sample Input

2
4
1 2 3 4
3
3 6 9

## Sample Output

Case #1: 2
Case #2: 3

hujie

## Source

BestCoder Round #38 (\$)


text:
# Question #8de97

Dec 10, 2016

That is not an identity.

#### Explanation:

Recall that

${\cot}^{2} x + 1 = {\csc}^{2} x$.

So, we can write

$\frac{1 - {\csc}^{2} x}{\csc} ^ 2 x = \frac{1 - \left({\cot}^{2} x + 1\right)}{\csc} ^ 2 x$

$= {\cot}^{2} \frac{x}{\csc} ^ 2 x$

Recall also that $\cot x = \cos \frac{x}{\sin} x$ and $\csc x = \frac{1}{\sin} x$.

This allows us to continue

$= \frac{{\cos}^{2} \frac{x}{\sin} ^ 2 x}{\frac{1}{\sin} ^ 2 x}$

$= {\cos}^{2} \frac{x}{\sin} ^ 2 x \cdot {\sin}^{2} \frac{x}{1}$

$= {\cos}^{2} x$

Which is not identically $\cos x$.

(${\cos}^{2} x = \cos x$ only when $\cos x = 1$ or $0$)

Dec 10, 2016

No. It is equal to ${\sin}^{2} x - 1$.

#### Explanation:

If we have $\frac{1 - x}{x}$, we can write it as $\frac{1}{x} - \frac{x}{x}$.

The same way, $\frac{1 - {\csc}^{2} x}{\csc} ^ 2$ can be written as $\frac{1}{\csc} ^ 2 x - \frac{{\csc}^{2} x}{{\csc}^{2} x}$.

This is equal to ${\sin}^{2} x - 1$.